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Abstract of the Dissertation 

Hydrodynamic Description of 
Dilepton Production 

by 

Kevin Dusling 

Doctor of Philosophy 

in 

Physics 

Stony Brook University 
2008 

The first part of this thesis focuses on the production of thermal 
dileptons from a hadronic gas at finite temperature. The rates 
are calculated by an expansion in Pion density and constrained by 
broken chiral symmetry and vacuum correlation functions, many 
of which have been measured by experiment. We focus on emission 
processes having two Pions in the final state. 

Next, follows a separate discussion on viscous hydrodynamics and 
its effect on p±_ spectra and elliptic flow. A non-central hydro- 
dynamic model of Au-Au collisions in 2+1 dimensions is simu- 
lated. Off-equilibrium corrections to the distribution can bring 
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about large changes in the differential elliptic flow, especially at 
higher p±. Also discussed is the shear viscous correction to dilep- 
ton production in a quark-gluon plasma (QGP) emanating from 
qq annihilation in the Born approximation. It is argued that 
a thermal description is reliable for invariant masses less than 
M max « {2tqTq) I '(rj I s) . Shear viscosity leads to qualitative dif- 
ferences in dilepton p± spectrum, which could be used to extract 
information on the thermalization time, viscosity to entropy ratio 
and possibly the thermalization mechanism in heavy-ion collisions. 

Finally, the dilepton rates used in this work are integrated over 
the space-time evolution of the collision region and compared to 
the recent results from the NA60 experiment at CERN and the 
PHENIX experiment at RHIC. The role played by chiral symmetry 
restoration in the hadronic phase and viscosity in the QGP phase 
is discussed. 
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Introduction 

1.1 Basics of Quantum Chromodynamics 

Quantum chromodynamics (QCD) is believed to be the theory which describes 
the strong interactions of quarks and gluons which are found in hadrons. The 
dynamics of both quarks and gluons are dictated by the QCD Lagrangian [1] 

1 f) 

L.QCD ~ --^^a + Y^3 e <^n^ a /3 

+ (i.i) 

/=! 

where we have used the following notation 

G% = d,A a v -d y Al + igt bc AlAl 

D$ = d^ + ig{A;T a f (1.2) 

1 



and A® is the gluon field having color index a and is a quark field having 
flavor index / and color index i. 

Unfortunately there is no calculational scheme which works well for all 
energies. For large momentum transfers asymptotic freedom states that the 
coupling constant becomes small, therefore allowing for a perturbative treat- 
ment. However, as the momentum transfer decreases the coupling becomes 
larger binding the quarks and gluons into hadrons. This leads to the property 
of confinement, whereby the force required to separate two quarks increases 
with their relative distance. 

Since perturbative calculations are only allowable at high energies due to 
asymptotic freedom non-perturbative methods have been developed in order 
to gain insight into QCD. One of the more well established non-perturbative 
approaches is lattice QCD. Only recently have first principal calculations by 
lattice QCD in the strong coupling regime become available, albeit with limi- 
tations. Since the metric in lattice QCD is Euclidean the calculation is limited 
to static properties. For example, it becomes very difficult, if not impossible, 
to calculate scattering amplitudes or transport coefficients. 

There have also been many non-perturbative methods developed based 
on effective theories. In order for these theories to represent nature they 
should contain the same symmetries of QCD. QCD with massless quarks and 
Nf flavors has an exact global flavor symmetry SUi{Nf) x SUpi{Nf) called 
chiral symmetry. This symmetry is spontaneously broken generating three 
(for SU(2) flavor) Goldstone bosons called the Pions. In the real world the 
quarks are massive and electromagnetism is present, so the flavor symmetry 
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is only approximate, leading to pseudo-Goldstone bosons having a small mass 
which can be calculated in the framework of chiral perturbation theory. This 
explicit breaking gives rise to the partially conserved axial current (PCAC) 
hypothesis, d^A't = f 7T m^ir a with A a ^ the axial-vector current. 

1.2 Heavy Ion Collisions 

Chiral symmetry which is spontaneously broken in the QCD vacuum is par- 
tially restored at finite temperature and/or density. In addition, as the tem- 
perature is increased from zero, it is thermodynamically favorable for there 
to be a phase transition from a resonance gas of hadronic bound states to a 
quark gluon plasma. One of the goals of the heavy ion collision program is to 
produce a quark-gluon plasma and study its properties. It has already been 
accepted by many in the heavy ion community that a quark gluon plasma 
(QGP) has been triggered at RHIC consisting of a strongly interacting, low 
viscosity fluid. 

In order to confirm these conclusions and quantify the properties of the 
QGP, a detailed study of heavy ion phenomenology is required. One of the 
main experimental observations that led to the conclusion of the low viscosity 
nature of the QGP is the large amount of collective flow of the produced 
particles and its interpretation as coming from a hydrodynamic expansion. 
Even though hydrodynamic behavior is able to explain a large amount of the 
available hadronic data it fails at a number of places, such as at high transverse 
momentum (jpx) and at forward rapidity. It is believed that the deviations from 
ideal hydrodynamic behavior could be explained by dissipative effects. 
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In order to quantify these assertions, viscous relativistic hydro dynamic sim- 
ulations have to be developed. The first order Navier-Stokes theory is plagued 
with difficulties (e.g. the parabolic nature of the equations permit acausal sig- 
nal propagation). In order to correct for this unsatisfactory behavior a number 
of second order theories have been developed. At this point in time there is 
still not a consensus in the heavy-ion community on which theory is appropri- 
ate to use. A full study of viscous hydrodynamics is not only imperative for 
making quantitative predictions on the properties of the matter produced at 
RHIC but also helps in our theoretical understanding of kinetic theory results. 
The importance of having a theoretical understanding of viscous relativistic 
hydrodynamics is not limited to heavy-ion collisions but is also necessary, for 
example, in cosmological simulations of the early universe. 

A second interesting phenomenological tool to study heavy ion collisions 
is electromagnetic probes. In contrast to hadronic observables which interact 
strongly throughout the entire evolution of the heavy ion collision, electromag- 
netic probes leave the medium without further interaction and therefore carry 
direct information on the time evolution of the system [2] . This is in contrast 
to hadronic observables which thermalize after the collision and thus provide 
information only on the late stages of the evolution. 

In theory, the electromagnetic spectral function of the quark-gluon plasma 
could be extracted from thermal photon and dilepton emission, which would 
in turn permit one to learn about its properties (e.g. transport coefficients, 
presence of bound states, etc.) In practice, however, this is not possible since 
the QGP dilepton yields are quenched by hadronic emission. Therefore, in 
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order to probe the QGP phase, there must be a solid theoretical understanding 
of the hadronic emission processes. 

There is a long history of experimental dilepton measurements [3] which 
we don't attempt to summarize here. In regard to dilepton measurements 
from heavy-ion collisions there were three experiments prior to the recent re- 
sults from NA60 and PHENIX. These three past experiments were the NA45, 
HELIOS-3 and NA38/50, which were all performed at the CERN SPS col- 
lider, and focused respectively on low, intermediate and high mass dileptons. 
All three experiments found an enhancement in the dilepton yields above ex- 
pected hadronic sources (which is comprised of a cocktail designed to describe 
the measured dilepton spectra in p-p and p-Be collisions). The quality of 
data however was limited. Dilepton measurements in general are much more 
difficult then measuring hadronic observables. Not only is there a large back- 
ground which must be rejected but the cross sections involved are also relatively 
small. These two facts together require high luminosity experiments in order 
to collect precision data. 

Let us discuss the low mass enhancement found at NA45. For central S-Au 
and Pb-Au collisions NA45 found an enhancement by as much as a factor of 
3-5 above known sources in the emission of di-electrons with invariant masses 
0.2 < M < 0.6 GeV. A number of theoretical explanations were given for this 
phenomenon including melting of the p due to chiral symmetry restoration 
[4, 5, 6, 7, 8]. The statistics were unfortunately too poor to confirm these pre- 
dictions. The upgrade from NA50 to NA60 consisted of a new vertex tracker, 
which now allows track matching in both coordinate and momentum space. 
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This leads to a considerable improvement in statistics and should allow one to 
discern between different theoretical approaches. 

Recently, the PHENIX experiment has also measured di-electron invariant 
mass spectra and found an enhancement by a factor as large as 7-8 above 
the cocktail for the most central collisions. At first glance it appears that 
this result may be inconsistent with the measurements by NA60. However, 
one must remember that the resultant yields must first be folded through the 
complicated detector acceptance which is specific to either PHENIX or NA60. 

In order to support these interesting experimental programs it is necessary 
to generate realistic dilepton predictions. In this direction we use a compre- 
hensive set of rates for dilepton production taking into account the symme- 
tries of QCD (e.g. broken chiral symmetry) at finite temperature and density 
integrated over the space-time history of relativistic viscous hydrodynamic 
simulations of the collision. 

1.3 Outline of this thesis 

This thesis is separated into a number of self contained parts. However, the 
last section on heavy ion phenomenology will rely on all the material presented 
throughout. 

The first part of this work focuses on the work done with my advisor, Ismail 
Zahed. In chapter two, we discuss the dilepton emission rates from a hadronic 
gas in thermal and chemical equilibrium. The rates take into account broken 
chiral symmetry in a consistent manner and rely on experimental data as 
input. The rates are treated in a density expansion and the effects of one and 
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two pions in the final state are explored. The new work consisted of evaluating 
the rates to second order in pion density, which include all hadronic processes 
involving two pions in the final state. This work is currently unpublished. 

Chapter three is a separate discussion on viscous hydrodynamic simula- 
tions. The work was done under the auspices of Derek Teaney and was pub- 
lished in [9]. We examine how shear viscosity changes the ideal hydrodynamic 
evolution and the effect it has on differential transverse momentum and elliptic 
flow spectra. 

Chapter four goes back to dilepton production, this time from a quark 
gluon plasma out of kinetic equilibrium. In collaboration with Shu Lin [10] 
we consider how shear viscosity modifies the leading order born qq dilepton 
production rates. 

In chapter five all the pieces are put together. The dilepton rates are 
integrated over the space time evolution presented in chapter three. Most of 
the work is done in kinetic equilibrium and was published in [11, 12, 13] and 
was done in collaboration with Ismail Zahed and with help from Derek Teaney 
regarding the hydrodynamic evolution and equation of state. A final section 
in chapter five discusses the role of shear viscosity on dilepton emission from 
both the QGP and hadronic phases. 
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Chapter 2 

"Master Equation" Approach to 
Dilepton Production 

2.1 Introduction 



It can be shown [14] that to lowest order in electromagnetic interactions and 
to all orders in strong coupling the differential rate for dilepton pair produc- 
tion can be expressed in terms of the correlation function of the hadronic 
electromagnetic current. 

When lepton mass is ignored the rate is given by 



dR Aa 2 1 
~¥q ~ 3(2vr) 3 ^ 



- 4 (q"q u ~ q 



(2.1) 



where 
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In the above equations q is the time-like four-momentum of the lepton pair, 
J^ m is the hadronic part of the electromagnetic current and (• • • )p stands for 
the thermal averaging at a temperature f3 = 1/T. 

In general, there are two ways in which the above thermal structure func- 
tion, Wfj, v (q), can be evaluated. The first is by kinetic theory. By inserting 
a complete set of states for each incoming component of the thermal density 
matrix the above equations can be shown to agree with relativistic kinetic the- 
ory reaction by reaction. This method of evaluation is not only cumbersome 
but also relies on many approximations, such as the choice of Lagrangian and 
coupling constants. 

A second approach, which is used here, is to relate the thermal structure 
function directly to spectral functions as was first done by Z. Huang [15]. From 
the spectral representation and symmetry the thermal structure function can 
be related to the absorptive part of the time ordered function 

W^( g ) = 1+ ^ /r Imn^( g ), (2.3) 

where 

np(?) = • J d 4 xe^(T {j; m (x)Jr(0)))e ■ (2-4) 

Let us now evaluate H e ™ not reaction by reaction as done in kinetic theory 
but instead in a low temperature expansion as was first done by Dey, Eletsky 
and Ioffe [16]. At low enough temperature the heat bath will be dominated by 
pions and therefore one can keep the first term in the expansion of the trace 
in the thermal averaging. We quote the result and leave the details for the 
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appendix. To leading order in (T/f n ) 2 one finds 

U^q, T) = np(g) - e (n£(?) - U A u (q)) , (2.5) 

where e = T 2 /6f% and U A ' V are the axial- vector and vector correlators. This 
is an example of how the dilepton production rate at finite temperature can 
be determined from measurable experimental data at zero temperature. The 
above result shows that the vector and axial-vector spectral densities mix at 
finite temperature. Chiral symmetry makes the statement that U v = IT 4 . To 
leading order in temperature this occurs when T = y/3f n « 160 MeV. 

Even though the above result was restricted to zero momentum pions it 
is general in the sense that it was derived from current algebra and PCAC 
alone. As the number of soft-pion fields emitted or absorbed grows the current 
algebra formulation becomes increasingly difficult. For this reason Weinberg 
[17] developed a method of calculating current algebra results using an effective 
Lagrangian formulation at tree level. By renormalizing the tree level results 
one could obtain corrections to the soft-pion theorems. A one loop calculation 
can still only describe data up to about 200 MeV above threshold. Two loop 
calculations are intractable since over 100 new low energy constants appear. 

A program that extends chiral symmetry consistently into the resonance 
region without the soft-pion restriction is discussed in the next section. 
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2.2 Chiral Reduction Formula 



The limitations of current algebra and chiral perturbation theory can be 
avoided by instead using an S matrix formalism. H. Yamagishi and I. Za- 
hed [18] have derived a single equation (coined the "Master Equation") that 
contains all of the low energy theorems of current algebra. 

The starting point for this program is an action I with its kinetic part 
invariant under local SUl(2) x SUr(2) that is gauged with external sources. 
Examples are two- flavor QCD or the nonlinear sigma model. 

For two-flavor QCD the action is given as 

I = J d A xq^ (id, + G,+ t£y + a'yTte) q 

Tfl I 

- — \ I d A xq (ml + s - i"/ 5 T a p a ) q 
mi J v 

~ Y*f d 4 xTi c (G, u Gn (2.6) 
where G is the gluon field strength tensor defined by 

G/_ IU = d^G v — d v G ^ — i [Gn, G v ] (2-7) 
By Noether's theorem currents are defined by 



J <*> - m (2 ' 8) 



where J = (V, A, f n a, f w ir) and = (v®, a a , s,p a ). The currents must 
satisfy the Noether's equations in the presence of sources (also known as the 
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Veltman-Bell [19] equations): 

VV^ + fl^ + Apf^ (2.9) 

+ a"V„ - U (ml + s) vr + Upa = (2.10) 

In the above equations = d^l + v^ is the vector covariant derivative and we 
have used the notation that A ac = e abc A b . Schwinger's quantum mechanical 
action principal 

5(/5 ou t|«in) = i(Pout\5I\a in ) (2.11) 

along with the completeness of asymptotic states leads to the Peierls-Dyson 
formula [20] 



w-^mry (2 - 12) 



The Veltman-Bell equations can now be recast into the following form 



X v + p— )S = (2.13) 
-bp) 

X A -(ml + s)j- p +p^jS = (2.14) 

(2.15) 
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where we have defined 



X. a v (x) = Vl c — — - + — - (2.16) 
X a A (x) = V a u c -^— + a a u c (x)-^— (2.17) 

It can be shown that Xy and are the generators of local SU(2) x SU(2). 

So far we have not considered whether chiral symmetry is present or ex- 
plicitly and/or spontaneously broken. The spontaneous breaking of chiral 
symmetry is expressed in terms of the following asymptotic condition on the 
axial- vector field 

A» - -/A<it,in(aO *° - (2-18) 

The above condition assumes the absence of any additional stable axial vector 
or pseudo-scalar resonances. For explicit chiral symmetry breaking we must 
also impose 



PAftx) - +/X< l t,i» *° - ±oc (2.19) 

In order simplify the incorporation of the above boundary conditions into 
the approach a modified action (I) and modified S-matrix (S) are introduced. 

I = I - ft J d'x (s{x) + l -a^x) ■ a M (a0) (2.20) 
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S = S exp 



ifl I d A x ( s(x) H — a^{x) ■ a M (x) 



(2.21) 



Also needed is a change of variables, p = J/ f w — V^a^. We now use = 
(u" a",s, J a ) as independent variables and J = (j^, jy, /^a", 7r) as modified 
current densities defined analogously as 



3(x) = 



5S 



(2.22) 



Under this new change of variables the Veltman-Bell equations 2.15 can be 
integrated upon introduction of a retarded and advanced Green function Qn t A 
yielding a relation between the pion field and the other currents. After some 
manipulation these relations can be written in the following form 



TTin, S 



d 4 yd 4 ze lk y(l + Kg R ) ac (y,z) 
1 



-iS (K7r in ) c (z) + iUJ c {z)S - 1 (VM, - J) c (z)-^ -\ 



■X c A (z)S 



(2.23) 



With the above master equation in hand it is easy to see a strategy in 
order to generate Ward identities. After reducing out a pion from a scattering 
amplitude the commutator on the LHS can be replaced by an operator of 
many functional derivatives acting on the scattering matrix on the RHS. These 
variations on the S matrix can be written as time ordered products by 



T*(0( Xl )---0(x n )) 



S" 



8<t>(x x ) ■ --S(f)(x r 



-S 



(2.24) 
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where O = -i&5S/54>. 



2.3 Leading Order Lepton Emission: I — > Firl + l' 

The rate of dilepton emission per unit four volume for particles in thermal 
equilibrium at a temperature T is related to the thermal expectation value of 
the electromagnetic current-current correlation function [21, 22]. For massless 
leptons with momenta pi and p2, the rate per unit invariant momentum q = 
Pi + T>2 is given by: 

dR — o? 1 T . . . 

IT = tt^ - 9 i — I o7F ImW (?) 2-25 

d 4 q 3n 3 q 2 1 + e<rl T 

where a = e 2 /4ii, T is the temperature and 

W F (q) = i J d^x e^Tr [ e - (H - F)/T T*J^(x) J M (0)] (2.26) 

where eJ M is the hadronic part of the electromagnetic current, H is the hadronic 
Hamiltonian and F is the free energy. The trace is over a complete set of 
hadron states. 

In order to take into account leptons with mass mi the right-hand side of 
Eq. 2.25 is multiplied by 



2m? . , 4m? , , , 
[1 + —±)(1 ±) 1 ' 2 (2.27) 



Even though there are various approaches to calculating production rates, 
they differ in the way in which the current- current correlation function in 
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Eq. 2.25 is approximated and evaluated. The approach taken here is to use a 
chiral reduction formalism in order to reduce the current-current correlation 
function in 2.26 into a number of vacuum correlation functions which can 
be constrained to experimental e + e~ annihilation, r-decay, two-photon fusion 
reaction, and pion radiative decay experimental data. 

For temperatures T < m n the trace in Eq. (2.26) can be expanded in pion 
states. Keeping terms up to first order in pion density yields [6] 

ImW F (g) = -3q 2 lmU v (q 2 ) + i / ^Wf (g, k) (2.28) 

with the phase space factor 



d 3 k 1 1 , 



The first term in 2.28 is the transverse part of the isovector correlator 
(0|T*VV|0) which can be determined experimentally from electroproduction 
data and gives a result analogous to the resonant gas model. At low and 
intermediate invariant mass the spectrum is dominated by the p(770 MeV) 
and p'( 14 50 MeV). 

The term linear in pion density (the second term in Eq. 2.28) can be related 
to experimentally measured quantities via the chiral reduction formulas [23]. 
It is shown in [6] that the dominant contribution comes solely from the part 
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involving two-point correlators which yields: 
WZ(q,k) = ^g 2 Imn y (g 2 ) 

J 7T 

- l(k + q) 2 lmU A ((k + q) 2 ) + (q^-q) 

J TV 
g 

+ ~Pi ((k ■ q? ~ mlq 2 ) Imn^g 2 ) x ReA fl (A; + q) + (q -> -q) 

(2.30) 

where ReA# is the real part of the retarded pion propagator given by 1/ (q 2 — 
m 2 + ie) and is the transverse part of the iso-axial correlator (0|T*j^j^|0). 
The spectral functions appearing in Eq. (2.30) can be related to both e + e~ 
annihilation as well as r-decay data as was compiled in [24]. 

It can be seen in Fig. 2.1 that the term linear in pion density decreases the 
rates from the resonance gas contribution for the mass region above the two 
pion threshold. However below the two pion threshold the only contribution 
to the rates come from the terms in Eq. 2.30. This is because the axial 
spectral density is integrated over all momentum in the thermal averaging 
(Eq. 2.28), which weakens the (k + q) 2 factor in Eq. 2.30 allowing the 1/q 2 
term in Eq. 2.25 to dominate at low q 2 . 
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0.2 0.4 0.6 0.8 1 1.2 1.4 

M [GeV] 



Figure 2.1: (Color online) The total integrated dimuon rates from a pion gas 
at T=150 MeV. The curve labeled "Res. Gas" shows the analogue of the 
resonance gas contribution (the first term in Eq. 2.28). The curves labeled 
"Vector" and "Axial- Vector" show the contributions from the respective spec- 
tral functions in equation 2.30. 
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2.4 Second Order Lepton Emission: / — > Frnvl^l 
2.4.1 Introduction 

We now keep terms up to second order in pion density 

ImW F (g) = -3g 2 Imn y (g 2 ) + -^ J rfvrWf(g, k) + ^- j dn^W^q, k u k 2 ) 

(2.31) 



with the phase space factor 



d 3 k 1 1 



The first two terms in the above density expansion were considered in 
the previous section. As long as the system is sufficiently dilute (i.e. k = 
n n /2m n f% <C 1) the expansion should converge rather quickly as long as no 
new thresholds open up. What we will find is that the 2ir contribution feeds 
into the low mass and low p± region where the zero and first order corrections 
do not contribute. This will also enhance the real photon rate (M 2 = 0) at 
small energy. 



2.4.2 Result 

In this section we quote the full on-shell Ward identity for W^. We note 
that we correct a number of typographic errors from the result quoted in [7]. 
We also discuss in detail which terms are kept in the numerical calculations 
and argue which terms can be safely neglected. First let us quote the terms 
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included in the analysis: 

J TT 

= 4 [9^u - {2kt + q) fl k lu ReA R {k 1 + q)} ImTf" (q, k 2 ) 

J IT 

+ (q -> ~q) + (h -fci) + (q : h -> -g, -fcx) (2.33) 



+ ^^(2fci + g) 1/ ReA i? (A; 1 + g^V^Im^^, fc 2 ) 

i 7T 

i [<T - + g) M (2fc 1 + q yReA R (h + g)] 
x e a3e e a3/ Im^{(A; 1 + g,fc 2 ) 

+ + qY(h + qYiZh + qf [ReA R (h + q)f 

J TV 

x e^e^lmBlUh + q^^ + ih^ -h) (2.34) 



Equation 2.33 contains the pion-spin averaged ir'y forward scattering am- 
plitude (iT^j). This quantity can be constrained from measured photon fusion 
data by crossing and is discussed in section 2.4.3. We have defined the term 
B in equation 2.34 as 

BfKh, h 2 )=zj d*xe ik ^(ir b ont (k 2 )\T* (jy*)jL(0)) !<(**)) (2.35) 

The pions in the above expression for B can be reduced out via the chiral 
reduction formula. Since most of the strength will come from the vector and 
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Figure 2.2: Schematic representation of the chiral reduction of B. The dashed 
line is a pion. 

axial- vector spectral densities we keep these terms only. 

J 7T 

+ (k 2 -> ~k 2 ) 

- ^ ef [g^kl - k llM k lv ] ImIT A (kf) (2.36) 

J 7T 

There are additional terms however which we should discuss. Most can be 
argued away by resonance saturation. One term which we should quote which 
could appear at higher mass is the following four point function 

j 2 k«k p 2 i J d i xd 4 z l d i z 2 e lk ^e^e- ik ^(0\T* (jL(^)j^(^i)j'V(x)ji(0)) |0) com , 

but we neglect it further in this analysis since we want to focus on the re- 
gion below the p mass. Experimental information about this term could be 
extracted from n — n scattering data [1] . 

For completeness we now quote the remaining terms of *W^ n . These are not 
included in our analysis since they can be argued to be small in the kinematic 
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regions which we are interested in. 



+ ^k?k{{2k x + qfReAnih + q)ReA R (k 1 )e a3e e^ e JmB a a $(k 1 , fc 2 ) 

+ (q -5) + {h -fci) + (g, fci -g, -A*) (2.37) 



+ l^(2fc 1 + g)^ 3e W^(^i,g) 

+ (5 -5) + (*i ~h) + (q, h -g, -fci) (2.38) 



(2Ai + qf [ImA*(g + A*)] 2 e a3e e o3/ ImT^ 6/ ((h + q), fc 2 ) 

+ (g--g) (2.39) 



l m £_x / rfW^^-^(7rl(fc 2 )|T*V 3 (x)V 3 ^(y)a(0)|<(fc 2 )) 



(2.40) 

-^g^k«e a3e ImiC™ a (q, k x + g, A*) + (A* - -k x \ 



+ Im-^(2fc 1 + qYih + g)"fc«ReA fl (A; 1 + q)e a3e iC™ a (q, fci + g, fci) 

J IT 

(2.41) 

- j^kfk% j d 4 xd 4 yd 4 ze ik ^ y - x) e iqy 

J n J 

x Imz« t (fc 2 )|T*j^(x)j^(y)V 3 (z)V 3 ^(0)K(A; 2 )) (2.42) 



2m 2 

— ^(2fc! + q) 2 ReA R (h + q) [ReA^fci + g) + ReA^A*)] 

J 71" 

x Im(7r o 6 ut (A; 2 )|a(0)|7rL(A; 2 )> 

+ (g -g) + (h -fci) + (g, fci -g, -fci) (2.43) 
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We now discuss why the above terms are neglected. First look at eq. 2.37. 
It is proportional to the principal value of the real part of the retarded pion 
propagator defined as, 



A R (k) = PP , 2 1 2 - i-KSgn(k°)5(k 2 - ml). (2.44) 



For on shell pions this term is proportional to PP^ = and therefore vanishes. 
Now look at eq. 2.38 where we have defined A as: 



A%l(h,q) ee J d*xd*ye^e^(7r b out (k 2 )\T* (j^(x)Vj(y)) ^ n (0)K(h)} 

(2.45) 

Making use of the chiral reduction formula one can reduce the incoming pion 
7rf n (0) with the result: 



[ImAnih + q)]" 1 lmA™(k u q) = ^ImB^, k 2 ) 

Jit 

+ -U e3 * (2h + q)^ k{ ReA^lmB^ih, k 2 ) 

Jit 

+ ^ [ ci 4 xe^Imz(7rl(A ;2 )|T*V3(a ; )V^(0)|^ 1 (A ;2 )) 



- -^{h + qf I d4xd4 v eigXeiklv 

JIT J 

x Imi(7rl(fc 2 )|T*Vj(x)j^( 2 /)j^(0)|7rf n (A :2 )} 
+ ^8 ae (2k 1 + q) a [ d*xe lklX 

JIT J 

x Im(7rl(A; 2 )|T*Vj(x)a(0)|7rf n (A; 2 )> (2.46) 
From reducing out the incoming pion we find that A oc ImA^(/ci + q). The 
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Figure 2.3: Diagrams of the matrix elements in eqns 2.40, 2.41 and 2.42 (left 
to right). The dashed line is a pion and the wavy line denotes a photon. 

term only contributes when the pion from the heat bath has the kinematics 
specified by the delta function. Due to the small amount of phase space this 
term will be suppressed compared to the terms in 2.33 and 2.34. The same 
argument can be made for neglecting eq. 2.39. 

The matrix elements appearing in terms 2.40, 2.41 and 2.42 are shown in 
figure 2.3 where we have defined the term C in eq. 2.41 as 



It turns out that these three contributions can be argued to be small. First 
note that eqn. 2.40 vanishes in the chiral limit. Furthermore, after reducing 
out the incoming pions the remaining vacuum spectral functions will mostly 
consist of (0|VW|0) and (0|AAcr|0) for which the resonance saturation is 
small. Equation 2.41 depends on C which after chiral reduction mostly reduces 
to the three correlators; (0|VAA|0), (0|VVi4|0) and (OlVVVlO) for which there 
is no s-channel cut through resonance saturation. Finally the matrix element 
in eq. 2.42 will mostly contribute in the four and six 7r range and higher. Since 
we are focusing our attention near the 2ir threshold it is safe to say that the 




x <7rl(fc 2 )|T*Vj 



(0)jUs)&.(y)l<(*a)> (2-47) 
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above processes can be neglected since the correlators contribute at higher 
mass. In addition these processes are Boltzmann suppressed in comparison to 
the resonance and one n final state reactions. 

The final expression, eq. 2.43, is a direct consequence of the way chiral 
symmetry is broken through (0|<r|0). The term can be related directly to 
the scalar form factor through Im(7r 6 (fc2)|(3"|7r b (/c2)) = ImFg(£ = 0) = and 
therefore vanishes. 

2.4.3 Pion Compton Scattering Amplitude 

We consider the reaction j(qi) + 7r (A;i) — > 7(52) + ^(£2) and define the 
Mandelstam variables to be 

s = (k 2 + q 2 ) 2 
t = (h-k 2 ) 2 

u = (k g -q 2 ) 2 (2.48) 

Let us express the total Compton scattering matrix element Ai as 

M = e 2 e*i( qi y 2 {q 2 )T, u (2.49) 

The pion compton scattering amplitude in the born approximation is by 
now a textbook example [25]. The three Feynman diagrams of figure 2.4 
contribute which evaluate for forward scattering {i.e. q\ = q 2 = q and ki = 
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9i 



32 gi 



<?2 gi 



k-2 ki 



hi k\ 



Figure 2.4: Tree level contribution to the n — 7 scattering amplitude. 



ki = k) to 



(2k + q) u (2k + q)^ , {2k - q) v {2k - q), 



rat 



u — mi 



29, 



/lis 



(2.50) 



where 



it 



(2.51) 



However, since we are always below threshold the amplitude is always real 
and does not contribute to the imaginary part of the amplitude. To go further 
we make use of the master formula approach of the 77? — > 771 reaction [18]. The 
pion compton scattering process for real photons (q 2 = 0) was examined in 
The dominant contributions from the chiral reduction of the ir / ~f* — > ir / ~f* 
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Figure 2.5: Schematic representation of the higher order correction from the 
chiral reduction of the 7 — n scattering amplitude. 

process is 



T^(q,k) = h2k + q),(-q 2 k l/ +(k-q)q u y 3e e eb3 ReA R (k + q))Tl v (q 2 ) 

+ (q -> -q) + (k -> -A;) + (q, k -> -g, -A;) 

+ ^e a3e e 6e3 (- M 2 + g M g,)n y (g 2 ) 

- ^e a3e e be3 (- 9lxv (k + qf + (k + g ) M (* + g)„) IT A ((A; + g) 2 ) 

+ (fc->-fc) 

(2.52) 



and corresponds to the diagrams in figure 2.5. Finally, we average over isospin 
(i.e. xl/35 ab ) obtaining 

lmTp(q,k) = ^-(2k + q) fl (-q 2 k u +(k-q)q u )ReA R (k + q))lmU v (q 2 ) 

+ (q — > -q) + (k — > -A;) + (g, k — > -g, -A;) 
4 

+ 7772 (^? 2 ~ ImIV(g 2 ) 

- ^ + g) 2 - (k + + q) v ) Imn A ((A; + g) 2 ) 

+ (k^-k) 

(2.53) 
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2.4.4 Soft Limit 

It is instructive to look at the above result for soft pions (k\, hi —> 0) in the 
chiral limit (m n = 0). To first order in pion density the result is given in [6] 



W^(g,0,0) = 12q 2 [lmIT y (g 2 ) - lulling 2 )] 



(2.54) 



To second order in pion density there are two terms which remain in the 
soft limit. They are the terms proportional to g^ u in equations 2.33 and 2.34. 
In this limit we have 



1 p 2 
— W^(q,0,0) = —g^ 

J it J n 



lmT<%(q, 0) + -e a3e e a3f lmB^{q, 0) 



j - 



4- -24 
2 



(lmIV(g 2 ) - ImIT A (g 2 )) (2.55) 



We can now substitute the above result into 2.31 and noting that each 
phase space integral leads to a factor of T 2 /24 we find 

ImIT em (g 2 , T) = ImII cm (g 2 , 0) - (e - L\ [lmll v (q 2 ) - ImII A (g 2 )] (2.56) 



where 



6/f 



(2.57) 



consistent with the result of Dey, Eletsky and Ioffe [16, 27]. 
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2.4.5 Dilepton Rates 

The phase space integrals in the rate equations were evaluated numerically. As 
discussed earlier we expect the dominate contributions to come from eqns. 2.33 
and 2.34 where 7^ 7 and B are given by equations 2.53 and 2.36 respectively. 

We now discuss the results which are shown in figure 2.7. The upper 
(lower) figure shows the virtual photon rates for q = (0.5) GeV. The red 
curve labeled "Res Gas" and the gray curve labeled "1st order" correspond to 
the first and second terms of equation 2.31. These are the results originally 
found in [6] and correspond to all processes with zero and one pion in the final 
state. At q = 0.5 GeV there is a huge enhancement in the low mass region, 
especially below the two pion threshold, and can be thought of as arising from 
processes of the type I — > Fit + e + e~ . However, this enhancement gradually 
disappears at lower values of q and is negligible at q = 0. 

Let us now discuss the new results of this work, which are contributions 
from processes of the type I — > Fnir + e + e~. The magenta curves labeled 
"By" and "-77iy" , show the contribution from the vector spectral densities in 
the expressions for B and T nl respectively. We should note that the ■jttv term 
is negative throughout and therefore its absolute value is shown. The blue 
curves labeled "— £>a" and "77Ta" show the corresponding contributions from 
the axial- vector spectral density. In this case u Ba" is negative throughout and 
its absolute value is shown. This first thing to note is that the contributions 
from Ba and ^tta,v are well below the first order contribution in all kinematic 
regions. Furthermore, there is even some cancellation between Ba and 771^ 
which makes these contributions even less significant. It should be pointed out 
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that these terms cannot be neglected in the soft pion limit if one wants to be 
consistent with current algebra as discussed in the previous section. However, 
at finite temperature they can be neglected. 

The dileptons emanating from the By term cannot be neglected however, 
and lead to a large enhancement below the two pion threshold at small mo- 
mentum where the first order term is suppressed. More precisely, at q = 
GeV this term makes up for all of the emission below the two pion threshold. 
By q = 1 GeV the second order term can be neglected since the first order 
term dominates. 




Figure 2.6: Example of one of the matrix elements contributing to B. 



In order to gain some physical insight into where the low mass and low q 
enhancement is coming from we show in figure 2.6 one of the main contributors 
to the terms in eq. 2.34. The density expansion blurs the relation between the 
spectral densities and kinetic theory so it is not possible to make a mapping 
from a given term in the xRF to the relevant physical reaction. Regardless, one 
can think of the imaginary part of the diagram as coming from the thermal 
decay of p,p',p"--- — > nhi, 7r&i, 7ra\ ■ ■ ■ — > 27r7*. The two final state Pions 
carry away most of the momentum leaving behind a low momentum photon. 
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Figure 2.7: (Color online) Differential virtual photon rates at q = and 0.5 
GeV. The curve labeled "Res. Gas" is the zeroth order contribution. The first 
order curve is the sum of the resonance gas and first order, , contributions. 
The remaining terms show the contributions from various parts of W^. By{A) 
and 7vry( J 4) shows the vector (axial- vector) contribution of B and T 17T . (Note 
that Ba and 77iy are negative thus absolute values are shown.) The solid 
(dashed) black lines shows the total rate evaluated at fi n = MeV (p^ = 100 
MeV). 



31 



Appendix: Hadronic Correlators at Finite Tem- 
perature 

Let us consider the vector spectral function at finite temperature when the 
heat bath is dominated by pions. Expanding the thermal trace one finds 

n^(g,T) = zjd'xe^TT(e-^ H - F ^T{v;(x)V^O)))^U^ u ( qi T = 0) 

+ J d'xe^ J jSkC B{h,T) ^ mmT W (a W' (0) ) l 7 ^)) 

(2.58) 

Let's start with the following matrix element: 

M% = (7r a (k)\T (v;(x)Vj(0)) \Ak')) (2.59) 
Using the LSZ reduction formula one out-going pion can be extracted yielding 
M% = -i(k 2 - ml) J d 4 ze lkz (Q\T (n a (z)V*(x)VJ (0)) |tt 6 (A;')) . (2.60) 
Using the PCAC hypothesis, 

= -4rd a A a a (z) , (2.61) 

the matrix element becomes 

M% = -i ik2 ~f* ] [ d*ze ik *(0\T (d«Al{z)V;{x)VM) lA*)) ■ (2.62) 

171 n J IT J 
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We will use the following form of the divergence of the time ordered prod- 
uct: 



d:(0\T(A a a (z)V;(x)VJ(0))\n m (k')) 
= (0\T(d a A a a (z)V;(x)VJ(0))\Ak')) 
+ S(z - Xo )(0\T ( [A a (z), V;(x)] ^'(0)) \Ak')) 
+ 6(z -0){0\T([A« (z),Vj(0)]v;(x))\Ak')) (2.63) 

Note that the term on the left hand size of the equality will vanish after 
integrating by parts, neglecting the surface term and finally taking the soft 
pion limit. The commutators in the above expression are simplified by the 
current algebra relations, which hold for equal times (xq = yo): 

[A a o(x)X(v)] = ie ahc A^{x)5\x-y) 

[A a o(x),4(y)] = ie ahc V<{x)5\x - y) 

[Vo(x)X(y)] = ie ahc V<{x)5\x - y) 

[V a (x),A b ^(y)] = ie abc Al{x)5\x-y) (2.64) 



Performing the integrals in the matrix element one is left with 

f % M% = e ail (0\T(A l ^x)Vj(0))\Ak'))e ikx 

= e^ l (0\T(V^x)Ai(0))\7T b (k')} (2.65) 

The incoming pion in the above matrix element can be reduced in a similar 
fashion. Once both pions are reduced the matrix element can be substituted 
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back into eq. 2.58 



n£,( g ,r) = n^(g,T = 0) 

d 3 k 



+ (/ 



(2vr) 3 2A;o 



n B (k /T) 



x ± [<0|T (^(^(O)) |0) - (0|T (^(x)K(O)) |0>] 

/ 7T 



The result for the integral in parenthesis is 

d 3 k ,. 



(2.66) 



r d s k T 2 



The final result can be written as 



U v (q, T) = (1 — e) n£(g, T = 0) + eUUq, T = 0) (2.6* 



where e = 

6/^ 



Appendix: Collisional Broadening of the p Meson 

In this appendix we briefly discuss how the p meson propagator is modified at 
finite temperature in the resonance saturation approximation. This procedure 
was carried out in [28, 29] but we include it in this appendix for completeness of 
this work. It will be interesting to compare how standard collisional broadening 
of the rho meson compares to the spectral function approach based on the 
chiral reduction formula discussed in full in the previous section. 
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The total width of the p meson is the decay plus collisional width 



tot J- decay 



+ r coll (2.69) 



The decay width is given by the imaginary part of the rho self-energy at 
one loop while the collisional width is given at the two loop level. Employing 
an effective Lagrangian for the p — tt interaction the decay width is given as 
[30] 

rde - = Jb ^ ~ Aml)m i 2n * (I) + x ] (2 - 70) 

An effective Lagrangian can also be used to calculate the two loop level 
p self energy. However, the various couplings can not always be measured 
precisely and instead we resort to resonance saturation in order to calculate 
the 7rp cross section. The thermal collision rate is 

/• OC rri 

r c oii = ^/ A( 8 XX) K i(vtyTX,00 (2-71) 

In the above equation 5^,57, is the spin/iso-spin degeneracy factor of the 
7T, p mesons. n p is the number density of p mesons and X(x,y,z) = x 2 — 
2x(y + z) + (y — z) 2 is the kinematical triangle function. We use the standard 
Breit-Wigner formula for the irp cross section 



7T v -\ „ „ B R T 2 



a 



iTF,F, _ „ (2.72) 



where 'R' refers to intermediate resonances (here we include 0, a\, ci2, to' , 
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s 1 --Sq - [GeV] 

Figure 2.8: The irp — > 7rp cross section. 

and 7r(1670) with properties taken from [31]) and q is the three-momentum of 
the p meson 

(2.73) 



Q = 7^ xl/2 ( s , m l> m l) 
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Chapter 3 



Viscous Hydrodynamics 



3.1 Introduction 
3.1.1 Motivation 

One of the first and most exciting observations from the Relativistic Heavy 
Ion Collider (RHIC) was the very strong elliptic flow in non-central collisions 
[32, 33]. The elliptic flow is quantified by the anisotropy of particle production 
with respect to the reaction plane v 2 

and can be measured as a function of pr, rapidity, centrality and particle type. 

The adopted interpretation of the v 2 measurements is that the medium 
responds as a fluid to the differences in pressure gradients in the x and y 
directions. The fluid then expands preferentially in the reaction plane and 
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establishes the observed momentum space anisotropy. This hydrodynamic in- 
terpretation is supported by the qualified success of ideal hydrodynamic mod- 
els in describing a large variety of data over a range of colliding systems and 
energies [34, 35, 36, 37, 38]. Nevertheless, the hydrodynamic interpretation of 
the flow results is not unassailable. A back of the envelope estimate of viscous 
corrections to hydrodynamic results [39] suggests that viscous corrections are 
actually rather large, i.e., the mean free path is comparable to the system 
size [40]. These estimates are best conveyed in terms of the shear viscosity 
to entropy ratio, rj/s. The conditions for partial equilibrium at RHIC are so 
unfavorable that at unless rj/s is small (say 0.5 or less), it is difficult to imagine 
that the medium would participate in a coordinated collective flow. 

From a theoretical perspective, it is difficult to reliably estimate rj/s in 
the vicinity of the QCD phase transition where the system is strongly cou- 
pled. Lattice QCD measurements of transport are hard (perhaps impossible 
[41, 42]) though recent efforts have lead to estimates which are not incompati- 
ble with the hydrodynamic interpretation of RHIC results [43, 44]. In a strict 
perturbative setting (where the quasi-particle picture is exact) rj/s is large 
~ l/g 4 . Nevertheless an extrapolation of weak coupling results to moderate 
coupling also leads to an t]/s which is perhaps reconcilable with the hydrody- 
namic interpretation [45, 46]. Finally, these perturbative estimates should be 
contrasted with J\f — 4 Super Yang Mills at strong coupling, where r]/s is 1/47T 
[47, 48]. Although J\f — 4 SYM is not QCD, the calculation was important 
because it showed that there is at least one theory where r]/s is sufficiently 
small that collective phenomena would be observed under conditions similar 
to those produced at RHIC. 
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From a phenomenological perspective one of the most compelling evidences 
for the hydro dynamic interpretation of RHIC flow results is the fact that the 
deviations from hydrodynamics are qualitatively reproduced by kinetic theory 
[49, 50]. In particular, kinetic theory calculations generically reproduce the 
flattening of ^(pt) at higher p?, and the reduction of elliptic flow at large 
impact parameters. Some aspects of these kinetic theory results can be under- 
stood by considering the first viscous corrections to the thermal distribution 
function [51]. These estimates motivated full viscous hydrodynamic simula- 
tions of the elliptic flow which will be performed in this work. Recently such 
viscous simulations have been performed by two other groups [52, 53]. A 
brief discussion of the history surrounding viscous relativistic hydrodynamics 
is given below. 

3.1.2 Viscous Hydrodynamics 

The Navier-Stokes equations describe viscous corrections to ideal fluid flow 
by keeping terms up to first order in gradients of ideal quantities [54]. The 
resulting equations are parabolic which permit acausal signal propagation [55]. 
For instance, the stress tensor instantaneously adjusts to any thermodynamic 
force, diUj. This is, of course, an unphysical picture since the stress tensor 
should relax to the thermodynamic forces over a typical collision timescale. 

One would therefore like a phenomenological theory that explains this re- 
laxation correctly. Much work has been done in this direction but there is 
still no completely satisfactory theory. Probably the most used model is that 
of Israel and Stewart [56, 57], but there are also others by Lindblom and Ge- 
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roch [58] , Pavon, Jou and Casas- Vasquez [59] and also by Ottinger and Grmela 
[60, 61] which is used in this work. In fact a wide class of models was developed 
by Lindblom and Geroch in two separate papers [58, 62] 

All of the above theories have the same behavior: they relax on small time 
scales to the first-order relativistic Navier-Stokes equations and have some 
generalized entropy which increases as a function of time. It was shown by 
Lindblom [63] that for a large class class of these second order theories, the 
physical fields should be indistinguishable from the simple Navier-Stokes form. 
To paraphrase Lindblom; any measurement of the stress energy tensor or par- 
ticle current on a time scale larger than the microscopic time scale will be 
indistinguishable from the Navier-Stokes theory. The differences between the 
causal theories and the acausal Navier-Stokes equations are indicative of the 
corrections quantitatively captured by the full kinetic theory. Nevertheless, 
the causal theories provide a qualitative guide to the magnitude of these cor- 
rections [64]. However, the form of these corrections implicitly assumes a good 
quasi-particle description which may not exist in a strongly coupled plasma 
[65]. 



3.2 The Hydro dynamic Model 

In the following section we outline the equations of motion for the hydrody- 
namical model used in the following simulations. We start by summarizing the 
well known first-order Navier-Stokes theory. Then we outline the equations re- 
quired for a second-order causal description of dissipative fluid dynamics. This 
is done assuming a boost invariant expansion as first proposed by Bjorken 
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[66], where the equations of motion are expressed in terms of the proper time 
r = y/t 2 — z 2 and the spatial rapidity r/ s = |ln|^. The cartesian coordinate z 
denotes the position along the beam axis while x,y label positions transverse 
to the beam axis. 

3.2.1 1 st Order Viscous Hydrodynamics - Navier Stokes 

Viscous hydrodynamics was originally formulated in the first-order Navier- 
Stokes approximation where the energy momentum tensor and baryon flux is 
a sum of their ideal and dissipative parts: 

T» v = cuV + (p + II) A"" + tt^ , (3.2) 
= nu» + j%, (3.3) 

where p, e, n and = (7, tv) are the pressure, energy density, baryon density 
and four- velocity of the fluid. We use the convention that 
giiv _ diag(— 1, +1, +1, +1) and therefore u^u^ = — 1. The dissipative terms, 
7r and jd depend on the definition of the local rest frame (LRF) of the fluid. A 
specific form of and can be found using the Landau-Lifshitz definition 
[54] of the LRF (u M 7r M ^ = 0), constraining the the entropy to increase with 
time and by working within the Navier-Stokes approximation (keeping terms 
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to first order in gradients only) resulting in 

tt^ = -n^u" + VV - - t^VpuP) , (3.4) 

o 

n = -cv^, (3.5) 

f d = (3-6) 

where k, t] and ( are the heat conduction, shear and bulk viscosities of the 
fluid with temperature T and chemical potential fi. The viscous tensor is 
constructed with the differential operator V M = A^ u d v where = g^ u + 
u^u u is the local three-frame projector, d^u" = d^u 1 " + T^u 7 is the covariant 
derivative and = l/2g l/a {d^g ai +d 1 g a ^—d a g lil ) are the Christoffel symbols. 

The transport coefficients in a quark-gluon plasma and also in the hadronic 
gas were studied in Refs. [39, 45, 46, 67]. It was found that the dominate 
dissipative mechanism was shear viscosity in both the QGP and hadronic 
gas. Bulk viscosity may however dominate in the transition region [68]. Heat 
transport can be ignored in the limit that yUg T which is the limit taken 
here. 

In the following work we will consider viscous effects in a quark-gluon 
plasma phase only. For this purpose we consider a constant shear to entropy 
ratio, rj/s = const and a massless gas p = l/3e. Future work will discuss 
viscosity in the mixed and hadronic phases. From this point on we will neglect 
the thermal conductivity. We keep the bulk viscosity in the equations for 
consistency, but always set ( = in any calculations. 
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3.2.2 2 order Viscous Hydrodynamics 

In order to render a second order theory it is necessary to introduce additional 
variables. These variables will relax on very short time scales to the standard 
thermodynamic quantities in the first order theory, but an evolution equation 
for them is still required in order to avoid acausal signal propagation. One such 
theory that has been used in a number of works was introduced by Israel and 
Stewart [56]. Instead we use a theory developed by [60, 61] due to its appealing 
structure when implemented numerically. However, as discussed above, all of 
these theories should agree (i.e., they all relax on short time scales to the same 
the first-order equations). 

We now summarize the evolution equations used in the current analysis 
following the mathematical structure outlined in Ref. [61]. We use a simplified 
version of the model for deviations of the stress energy tensor close to equilib- 
rium. The new dynamical variable that is introduced is the tensor variable 
which will later be shown to be closely related to the velocity gradient tensor, 
Ti^ u . The tensor variable c^ u is conveniently defined to have the property 

c^ u u v = u M , (3.7) 

and the energy momentum tensor is given by 

= (e - u ot P a %)u'V + . (3.8) 

The explicit form of the stress tensor P^ is given in [61] and has a fairly 
complicated form. The discussion is simplified by considering small deviations 
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from local thermal equilibrium and working in the local rest frame so that the 
stress tensor can be approximated as 

T l L j RF = p(S ij - a#) , (3.9) 

where a is a small parameter related to the relaxation time. The equations 
of motion are dictated by conservation of energy and momentum which is 
given by d^T^ u = 0. In addition, an evolution equation for the generalized 
mechanical force tensor is needed and is given by [61] 

u x (d x c^ v - d^cxu - d u c^x) = —c^u - — , (3.10) 

where c and c are defined as the isotropic and traceless parts of the tensor 
variable c^ u defined as 

Cf,u = ^(c^-l){r] fiu + u fl u u ) , (3.11) 
c i^iu 4~ UfjUy = c^ u + c^ u . (3.12) 

In the limit that the relaxation times (r ,T2) are very small the evolution 
equation yields 

c ij = r 2 (<V + d jU * - p ij d k u k ) + 2 -r^d k u k . (3.13) 

Substituting the above equation into T^ RF and comparing the result to the 
Navier-Stokes equation (3.6) the bulk and shear viscosities can be identified 
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as 

7] = r 2 pa , 

C=\r Q pa. (3.14) 

In the model proposed by Ottinger [61] the quantity a is related to the equa- 
tion of state, but in the linearized version it is simply treated as a constant 
parameter related to the relaxation time. We fix a = 0.7 in all calculations, 
which then fixes the relaxation times (t2, t ) as a function of rj and (. 

It is natural to ask what is the effect of the relaxation time on the theory. 
In some sense this was already answered by Lindblom [63]. He showed that the 
physical fluid must relax to a state that is indistinguishable from the Navier- 
Stokes form. Therefore we expect the physical velocity gradients to agree with 
those given by the auxiliary tensor variable &~ v as in Eq. (3.13) . We expect 
higher order gradient terms to be necessary when there are large deviations 
between any observable computed using the physical fields or the auxiliary 
field c^ u . This will be used as a gauge in order to find the limit of applicability 
of any hydrodynamic calculations. 

1+1 Dimensions 

We now outline the equations of motion for the stress-energy tensor and 
the generalized mechanical force tensor assuming a boost-invariant expan- 
sion as well as azimuthal symmetry with arbitrary transverse expansion. It 
is easiest to work in polar coordinates (r, r, <j), rj) and since there is no de- 
pendence on <fi or rj the four-velocity can be expressed as = (7,71V, 0,0) 
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where 7 = J- In this coordinate system the metric tensor is given by 

<T = diag(-l,l,l/r 2 ,l/r 2 ). 

The first two equations of motion are given by the conservation of energy 
and momentum, d^T^ = for v = r and v — r. (Due to boost invariance 
and azimuthal symmetry the v = rj and v = equations are trivial.) 



d T T 00 + d r T m = — (T 00 + P 33 ) - -T 01 (3.15) 
r r 

d T T 01 + d r T n = — T 01 - i(T n - P 22 ) (3.16) 
r r 



where P 22 = r 2 P 22 and P 33 = r 2 P 33 . The evolution equations for the generic 
mechanical force tensor are: 



9 r c n + vd r c u - -[(1 - c 11 )^ 1 + c 01 <9 rM °] = — c 11 - — c 11 
7 7 r o 7 r 2 

«9 T £ 22 + W 9 r c 22 + -(£ 22 - c 11 ) + -c 10 = — t 2 - —d 22 
r r 7r 7r 2 

s33 1 „.q -33 1 2,_ 33 0\ %V 1q — 1 ^33 1 £ 33 



<9 T c 33 + ?;«9 r c 33 + -(c 33 + c uu ) - — c lu = — c - — c' 
r r 7T 7r 2 



(3.17) 



where c 22 = r 2 c 22 and c 33 = r 2 c 33 . 



1+2 Dimensions 

We now consider the 1+2 dimensional case without azimuthal symmetry but 
still having longitudinal boost invariance and use a coordinate system whereby 
the coordinates transverse to the beam axis are cartesian, (r,x,y,r)). Since 
there is no dependence on t] the four-velocity can be expressed as = 
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7(1, v x , v v , 0) where 7 = , = = . In this coordinate system the metric tensor 
is given by g^ u = diag(— 1, 1, 1, 1/r 2 ). 

In this coordinate system the first three equations of motion are given by 
the v — t, x, and y components of the conservation law d^T^ = 0: 



d T T 00 + d x T 01 + d y T 02 = ^(T 00 + r 2 P 33 ) (3.18) 
d T T 10 + d x T n + d y T 12 = — T 10 (3.19) 
d T T 20 + d x T 21 + d y T 22 = —T 20 (3.20) 



The evolution equations for the generalized mechanical force tensor are: 



(9 T + + ^<9 y )c 22 + 2[(c 22 - 1)^ + c 21 <9, 

2 

(3 T + v x d x + v y dy)c m + -(c 33 - 
(d T + + Vydy)c 12 + c 12 (9^ + <9 y ^) + (c 22 - 1)9^ + (c 11 - 
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(3.21) 



Initial Conditions 



The hydrodynamic simulation is a 2 + 1 dimensional boost invariant model 
with an ideal gas equation of state p = -g e. The temperature is related to the 
energy density with the Nf = 3 ideal QGP equation of state. We have chosen 
this extreme equation of state because the resulting radial and elliptic flow are 
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too large relative to data on light hadron production. Thus, this equation of 
state will estimate the largest elliptic flow possible for a given shear viscosity. 
We note that for any non-central collision we have chosen a default impact 
parameter of b=6.5 fm. 

Aside from the equation of state, the hydrodynamic model is based upon 
reference [35]. At an initial time Tq = 1 fm/c, the entropy is distributed in 
the transverse plane according to the distribution of participants for a Au-Au 
collision. Then one parameter, C s , is adjusted to set the initial temperature 
and total particle yield. Specifically the initial entropy density in the transverse 
plane is 

s(x,y,r ) = ^^-, (3.22) 
t ax ay 

where is the number of participants per unit area. The value C s = 15 
closely corresponds to the results of full hydrodynamic simulations [35, 36, 37] 
and corresponds to a maximum initial temperature of To = 420 MeV at impact 
parameter b = 0. With the entropy density specified the energy density can 
be determined. This requires inverting the equation of state. 

In a viscous formulation we must also specify the viscous fields, i.e. the 
& v in the second order setup. Following the general philosophy outlined in 
Section 3.1.2 we will choose the d a> such that the stress tensor deviations are 

V, = -7i (V M u„) n = -cvx = ( 3 -23) 
Since at time r Q the transverse flow velocity and the longitudinal flow velocity 
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is Bjorken this means that at mid rapidity 



1 2 

n xx = iTyy = - - n zz = - rj d z u z IT = (3.24) 



To achieve this condition we first rewrite the flow equations for small 
and vanishing transverse flow. The Cjj equations become 

d T c 11 = — , (3.25) 

TO T 2 

^22 °22 

d T c 22 = — , (3.26) 
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(3.27) 



In writing this we have used the fact that for small velocity c 00 ~ — u°u° . 
Then looking for the quasi stationary state we set the time derivatives to zero, 
and use the relations c 1 ^ = \c\ 5 l i and c u = c u + <? J to find that 



c" = \^-\^, (3.28) 
3 r 3 r 



2r 2r 2 



c = 37-37- < 329 > 



33 2 to , 4r 2 



c 33 = ■=- + — ■ (3.30) 
or or 



3.3 Hydro dynamic Results 

The equations outlined in the previous two sections were integrated numeri- 
cally using the initial conditions described above. In this section we now show 
the results of the simulation. Before showing the results of the 2+1 dimen- 
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sional simulation we outline some of the main physics points using results from 
the 1+1 dimensional case. 

Fig. 3.1 shows the energy density per unit rapidity (upper) and the trans- 
verse velocity (bottom) at various times for both ideal hydrodynamics and for 
finite viscosity (rj/s = 0.2). The effect of viscosity is twofold. The longitudinal 
pressure is initially reduced and the viscous case does less longitudinal pdV 
work as in the simple Bjorken expansion [39]. This means that at early times 
the energy per rapidity decreases more slowly in the viscous case. The reduc- 
tion of longitudinal pressure is accompanied by a larger transverse pressure. 
This causes the transverse velocity to grow more rapidly. The larger transverse 
velocity causes the energy density to deplete faster at later times in the viscous 
case. The net result is that a finite viscosity (even as large as rj/s = 0.2) does 
not integrate to give major deviations from the ideal equations of motion. A 
preliminary account of this effect was given long ago [69]. 

We now present results of the 2+1 dimensional boost invariant hydrody- 
namic model. Fig. 3.2 shows contour plots of the energy density per unit 
rapidity in the transverse plane at proper times of r = 1,3,6,9 fm/c. The 
initial conditions (r = 1) is taken from the Glauber model discussed before. 

Fig. 3.3 shows contour plots of the transverse velocity at the same times 
of r =1, 3, 6, 9 fm/c. At r = 1 the figure is blank since the velocity in 
the transverse plane is zero as set by the initial conditions. By looking at 
the contours of constant v/c one can see that a finite viscosity increases the 
transverse velocity. 

Since we are interested in elliptic flow which originates from the initial 
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Figure 3.1: (Color online) Plot of the energy density per unit rapidity (top) 
and of the transverse velocity (bottom) at times of r = 1,3,6,9 fm/c, for 
rj/s = 0.2 (solid red line) and for ideal hydrodynamics (dotted blue line). 
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Figure 3.2: (Color online) Contour plot of energy density per unit rapidity in 
the transverse plane. The contour values working outward are for r = 1 fm/c: 
15, 10, 5, 1, 0.1, for r = 3 fm/c: 10, 5, 1, 0.1, for r = 6 fm/c: 3, 2, 1, 0.5 and 
for t = 9 fm/c: 0.5, 0.375, 0.25, in units of GeV/fm 2 . 
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Figure 3.3: (Color online) Contour plot of transverse velocity, v± = y/v% + v%. 
The inner most contour is for v± = 0.1 and increases in steps of Av± = 0.15. 
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spatial anisotropy of the collision region it is useful to see how the spatial and 
momentum anisotropy develop in time. We therefore look at the following 
three quantities [70]: 



((y 2 + x 2 )) 

//rpxx _ rpyy\\ 
l(j<xx _|_ JtyyU 



<M> = — (f^) — ( 3 - 31 ) 

where the double angular bracket ((•••)) denote an energy density weighted 
average. The spatial ellipticity (e x ) is a measure of the spatial anisotropy 
as a function of time. The spatial anisotropy is what drives the momentum 
anisotropy (e p ). This quantity can be thought of as characterizing the p 2 -, 
weighted integrated elliptic flow [71]. The final quantity ((vt)) is the average 
radial flow velocity. All three of these quantities are plotted in fig. 3.4 for 
r)/s=Q.2, 0.05 and 10" 6 . 

As already shown in the 1+1 dimensional case the finite viscosity case 
does less longitudinal work. The longitudinal pressure is reduced while the 
transverse pressure is uniformly increased in the radial direction, i.e. gives 
no addition v-i component. This causes the transverse velocity (as seen in 
{{vt}) and fig. 3.3) to grow more rapidly while e p lags behind the ideal case. 
Furthermore, the larger radial symmetric transverse velocity causes a faster 
decrease in the spatial anisotropy. This further frustrates the build-up of 
the momentum anisotropy e p . We therefore expect to see a decrease in the 
integrated v 2 as the viscosity is increased. This is indeed the case as will be 
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shown. However, this effect is small compared to the change in v% from use of 
the off-equilibrium distribution function. 

3.4 Freezeout 

As discussed in the introduction, ideal hydrodynamics is applicable when 
^mfp *C L where L denotes the typical system size. When dissipative cor- 
rections are included, one must remember that the Navier Stokes equations 
are derived assuming that the relaxation time tr is much smaller than the 
inverse expansion rate, trO^u^ <C 1. Therefore, in the simulations we deter- 
mine the freezeout surface by monitoring the expansion rate relative to the 
relaxation time using a generalization of the freezeout criteria first proposed 
in [72, 73] and later in [74]. 

Specifically, freezeout is signaled when 1 

V -dX ~ \ ■ (3-32) 

This combination of parameters can be motivated from the kinetic theory 
estimates [75]. The pressure is p ~ e (v% h ) with (u| h ) the typical quasi-particle 
velocity and e the energy density. The viscosity is of order rj ~ e (v^ h ) tr with 
tr the relaxation time. Thus the freezeout condition is simply 

-dX ~ T R d^ ~ ~ . (3.33) 



In actual simulations we take (r]/p) d^u 11 = 0.6 for most runs (see below). 
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Figure 3.4: (Color online) Time evolution of the spatial ellipticity e x , the 
momentum anisotropy e p , and the energy density weighted transverse flow 
(i v ±))i see Eq. 3.31. 

56 



In the model we are considering r\jp = ctT2 with a = 0.7 as described in 
Section 3.2.2. 

The value of | can be considered as a parameter chosen to be smaller than 
one. The point is that as the above quantity becomes large the Navier Stokes 
approximation is no longer applicable and the simulation should freezeout. At 
this point one would need to include further higher order corrections in the 
gradients or switch to a kinetic approach. 

It is also convenient to have a definition for an analogous freezeout surface 
in the case of ideal hydrodynamics. One can think of keeping the freezeout 
surface fixed as rj/s is taken to zero. Dividing the freezeout criterion by f]/s 
and using s = (e + p)/T ~ 4p/T we define 

X=^d,u\ (3.34) 

which involves only quantities in the ideal simulation. This is a separate 
freezeout parameter independent of the viscosity. We should point out that the 
ideal freezeout conditions becomes more complicated in a hadronic resonance 
gas phase. For example, the simple temperature dependence in x is modified 
due to the rho resonance peak in the it — tt cross section [76]. This will be 
considered in more detail in a future work. 

We show in fig. 3.5 contour plots of the freezeout surface for fixed x from 
both ideal (upper plot) and viscous hydrodynamics (lower plot). For fixed 
X the freezeout surfaces remain approximately the same in both cases. The 
freezeout surface from now on will be specified by x i n order to facilitate a 
comparison between the ideal and viscous cases when comparing spectra. 
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We have typically chosen x an d r]/s in order that ^d^u 11 = 0.6. Thus in 
Table 3.1 for rj/s — 0.2 we have x = 3.0 and = 0.6. However, for 

rj/s — 0.05 the freezeout parameter is x = 12 giving an unphysically large 
surface. This would normally not be the case in a more realistic model with 
a phase transition present, since in the hadronic phase the viscosity goes like 
7) ~ The change in scaling with temperature would cause the system to 
freezeout soon after hadronization. We plan on quantifying this statement in 
a future work. We therefore use (r]/p)d fl u tJ ' = 0.225 when t]/s = 0.05 giving 
X = 4.5. The thin solid curve in the lower plot of fig. 3.5 shows this particular 
freezeout contour. In table 3.1 we summarize the freezeout parameters used 
throughout this work. For a given r]/s the most physical choice of freezeout 
parameter x is selected such that (t]/p)d ll u ,J ' « 0.6. However, if the viscosity 
becomes so small that the volume becomes unphysically large (such as for 
rj/s — 0.05) we set x = 4.5 as a maximum. These three physically motivated 
parameter sets are given in bold in the table. 

We should stress that the freezeout surface taken in this work is different 
from the typical constant temperature surface used in many hydrodynamic 
simulations. From fig. 3.5, one can see from the temperature map that the 
surface is not an isotherm and actually spans a very wide range of tempera- 
tures. The freezeout surface is understood by examining the expansion rate in 
Bjorken geometry 

dX = d T u T + — + d x u x + d y u v . (3.35) 
The resulting surface is due to a competition between the first two terms in 
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3.35 at early times and the last two terms at later times. 



Tj/s 


29„u" 

p A 1 


X 


0.05 


0.6 


12.0 


0.05 


0.225 


4.5 


0.05 


0.15 


3.0 


0.2 


0.9 


4.5 


0.2 


0.6 


3.0 


0.133 


0.6 


4.5 



Table 3.1: Freezeout parameters used throughout this work. For a given 
rj/s the most physical choice of freezeout parameter \ is selected such that 
(r]/p)d IM u >1 ~ 0.6. However, if the viscosity becomes so small (such as for 
rj/s — 0.05) that the volume becomes unphysically large (see text for dis- 
cussion) we set x = 4.5 as a maximum. These three physically motivated 
parameter sets are in bold. 



3.5 Spectra 
3.5.1 Anisotropy 

Before computing the differential spectrum we will compute the momentum 
anisotropy as a function of time. The momentum anisotropy A 2 (which differs 
from v 2 by the placement of averages) is defined as 

A ^ = (pI) ~ (P 2 y ) _ S n - S 22 ^ 
(Pi) + ipy) S n + S 22 ' 

where S l i is the sphericity tensor and can be related to the hydrodynamics 
fields (i.e. ,71"^, II) and moments of the ideal particle distribution function. 
From a theoretical perspective, A 2 is preferred because it is almost independent 
of the details of the particle content of the theory [71]. 
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Figure 3.5: (Color online) Contour plot of various freezeout surfaces for central 
Au-Au collisions. Top: Surfaces from ideal hydrodynamics where the freezeout 
condition is set by the parameter x=1.5, 3 and 4.5. Bottom: Corresponding 
viscous solution where rj/s was fixed by the condition ^d^u^ = 0.6. The thin 
solid black curve shows the contour set by -<9 M w M = 0.225 for comparison. 
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We plot A 2 in the following manner. At a given proper time we integrate 
over the surface of constant Xi which has developed by time r. The remaining 
part of the surface is fixed by integrating over the matter which has not frozen 
out (x < Xf.o.) a t fixed proper time. This can be thought of as a freezeout 
surface with a flat top at time r. As time moves forward eventually all of the 
matter is frozen out over a surface set by constant \ yielding a constant A 2 . 

Figure 3.6 shows A 2 for four different freezeout surfaces. The figure on 
the top shows the results using only the ideal contribution to the sphericity 
(regardless of if viscosity is present). This will be analogous to using only the 
ideal particle distribution function when generating the spectrum. First look 
at the solid black curves which are generated using ideal hydrodynamics and 
a specified %■ F° r a larger value of % a larger space-time region is evolved 
by hydrodynamics producing a larger elliptic flow or A 2 . The true ideal case 
where hydrodynamics is universally applicable is given by \ = 00 • We see that 
for x — 4.5 most of the elliptic flow is reproduced. 

In order to assess the role of viscosity we first look at the figure on the 
top. The dashed curves show A 2 for 77/s = 0.05 and 77/s = 0.2 without 
including viscous corrections to the distribution function. (For clarity, these 
curves are shown only for \ = 3.0 and x = 4.5.) Without the corrections to 
the distribution function the viscous corrections to A 2 are modest. The lower 
figure shows the analogous plot, this time including the viscous corrections to 
the distribution function. The corrections are much larger and we therefore 
expect the viscosity to decrease the integrated elliptic flow. 
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Figure 3.6: (Color online) A<i (defined in Eq. 3.36) as a function of r. The solid 
black lines show the ideal result for \ =1.5, 3.0, 4.5 and 00. Also shown in the 
bottom and top figures respectively are the viscous results with and without 
including the viscous correction to the distribution function, for % = 3-0 and 
4.5 and r]/s — 0.2 (dashed green curve) and for 77/5=0. 05 (dotted blue curve). 
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3.5.2 Spectra 

The thermal pr and differential v 2 spectra of particles are generated using the 
Cooper-Frye formula [77] given by 

The thermal distribution function used in the Cooper-Frye formula above also 
needs to include corrections due to finite viscosity. We therefore write / = 
fo+Sf where f Q is the ideal particle distribution and 5f is the viscous correction 
given by 



7i> + -nA 



(3.38) 



5 

For boltzmann statistics f (l + f ) is replaced by f . The elliptic flow is defined 
as the weighted average of the yields with cos(20): 

dN 



r ds- dN 

J—n t dypxdpTC 



v 2 (p T ) = (cos(20)) PT = J ~" f ; ligp^g , (3.39) 



where <fi is the angle between the decaying particle's momentum (pt) and the 
azimuthal angle of the collision region. 

A typical freezeout surface for x = 3 at an impact parameter b=6.5 is 
shown in fig. 3.7. Color gradients show the temperature profile on the freezeout 
surface and as noted before the surface is not necessarily an isotherm. 

Differential pt spectra for massless particles are shown in fig. 3.8 for two 
different freezeout surfaces: x = 3-0 (top) and x — 4.5 (bottom). In both plots 
the ideal case is shown by the solid red line. First we discuss changes to the 
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Figure 3.7: (Color online) Freezeout surface for semi-central (b=6.5) Au-Au 
collisions for rj/s — 0.2 and \ = 3-0. 
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spectra brought about by modifications to the equations of motion by looking 
at the spectra generated with the ideal particle distribution (f Q only). For 
both values of viscosity and both freezeout choices a hardening of the spectra 
is observed. This is expected since viscosity tends to increase the transverse 
velocity. 

The effect from the viscous corrections to the distribution function are 
more subtle. At earlier times the transverse flow has not fully developed and 
the longitudinal pressure is reduced while the transverse pressure is increased 
[51]. This is a consequence of the fact that the shear tensor is traceless. The 
increase in transverse pressure leads to a hardening of the spectrum after 
integration over the space-time freezeout surface. This is the case for \ = 3 
even though the corrections are small. At later times the larger transverse flow 
alleviates some of the longitudinal shear. When the hydro is finally in a full 
3D expansion, the viscous correction tends to reduce the transverse pressure. 
This changes the sign of the viscous correction term. This is seen for \ — 4.5 
where the viscous corrections soften the spectrum slightly. 

As discussed above, any observable created by using the auxiliary variable 
& v should agree with the results using the physical velocity fields. Therefore 
we also show the viscous corrections calculated using the physical gradients 
(denoted by Sfc), i-e., in the local rest frame the 7r u is approximated by 

7r« = -7]{d l u j + <9V - -5 ij d lU l ) , (3.40) 

3 

when computing 5f. 

Overall, the corrections to the spectra are small so it is hard to see any 
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differences between the two calculations. This will not be the case for the 
differential elliptic flow where this comparison will be more important. 

Figure 3.9 shows the differential elliptic flow using the same parameter set 
from the px spectrum. The solid red curves shows the ideal spectrum and, as 
expected, a larger elliptic flow is generated for x = 4.5 compared to x = 3 since 
a larger fraction of the space-time volume is described by hydrodynamics. 

The viscous correction to the equations of motion causes only a small 
change in the elliptic flow as seen by comparing the results at finite viscosity 
using f only with the ideal case. For \ — 3 the change is almost negligible. 
For x — 4.5 deviations are on the order of 2% at px = 2 GeV. 

Including the viscous corrections to the distribution function can bring 
about large changes in the elliptic flow. We show the corrections due to the 
auxiliary variable by 5f n and those from the gradients by 8 fa and we expect 
the two results to agree. When the two results start to diverge the gradient 
expansion is no longer valid and a kinetic description is really required. 

Based on our discussion in section 3.4 the viscosity is what sets the freezeout 
surface. For rj/s — 0.2 we find that x = 3 (upper figure). In this case the 
viscous corrections are large but can only be trusted up to px ~ 1 GeV. We 
also show for comparison the spectra for rj/s = 0.05 which can be trusted past 
2 GeV. For r]/s — 0.05 we take x = 4.5 for reasons discussed in section 3.4. 
Again, the viscous correction decreases the elliptic flow as a function of px- 
Also shown are the spectra for rj/s = 0.2 and the corrections are larger. In 
both cases the spectra can be trusted past px = 2 GeV. 



66 




0.01 



0.5 



1.5 



Pt (GeV) 



Figure 3.8: (Color online) Differential transverse momentum spectra for Au- 
Au collisions at b=6.5 fm. The upper plot is for freeze-out parameter \ = 3 
and the bottom for x = 4.5. In both plots the ideal case is shown by the 
solid red curve. Then the viscous case is shown without including the viscous 
corrections to the distribution function and is denoted by f . The addition of 
the viscous correction to the distribution function is generated in two different 
ways. 5f n is calculated using the auxiliary tensor while 5 fa is calculated 
using the physical gradients i.e., = —r]{d IJ 'd v ). 
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Figure 3.9: (Color online) Differential v% spectra for Au-Au collisions at b=6.5 
fm. The resulting curves are generated in the same way as described for the 
Pt spectra in fig. 3.8 
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3.6 Conclusions 



In summary we now make several conclusions regarding the effects of shear 
viscosity on heavy ion collisions. 

We first recall the setup. The paper is restricted to an ideal gas equation 
of state p = |e and sets the initial non-equilibrium fields to the value expected 
from the navier stokes equations 7r u = —rj {d l u^). The initial distribution 
of entropy density follows the distribution of participants. (This could be 
changed to a Color Glass Condensate model initial conditions [78].) The paper 
simulates a fluid model based on [61] which is similar but differs from that 
of Israel and Stewart. However all models should ultimately agree on the 
magnitude of viscous corrections provided the viscosity is sufficiently small. 

Several technical notes warrant discussion here. An algorithm for a reliable 
solution of the viscous model was developed by Pareschi [79]. For small enough 
relaxation times the auxiliary fields 7r u should relax to the form expected from 
the Navier-Stokes equation 7r y ~ —^(d 1 ^). Generically, relaxation models 
for viscosity have long time parameters (the shear viscosity r/ in this case) 
and short time parameters. In the model considered here, a is the short 
time parameter while in the Israel-Stewart theory this short time parameter 
is Tj/[(e-\-p)r^\. These short time parameters can be constrained by the /-sum 
rule [80, 81, 82]. In general the results should not depend on these short time 
parameters. 

We now summarize our physical results. The integrated viscous corrections 
to the flow are small. This was seen in both the hydrodynamic fields and also 
in the differential and integrated elliptic flow when the thermal distribution 



69 



function was restricted to the ideal form. (The remainder of this paragraph 
discusses only results with this restriction.) For the integrated Vi this is seen 
in the upper plot of fig. 3.6 where A2 is shown for ideal runs and viscous runs 
at r]/s = 0.05 and 0.2. Corrections due to the modified flow pattern are also 
small in the differential v 2 spectrum as seen in fig. 3.9 by comparing the ideal 
and viscous runs (again with f only.) Although there is the possibility for 
the elliptic flow to be modified from variations in the freezeout surface across 
different runs this was minimized by freezing out on contours of constant %■ 
One can see from fig. 3.5 that the space-time freezeout contours are about 
the same at zero and finite viscosity. The fact that only small changes in the 
fields are seen when including viscosity is not surprising. The time scale of 
any heavy ion collision is much shorter than the time needed for dissipative 
effects to integrate and become large. 

Even though viscosity does not modify the flow strongly we have shown that 
there are still large corrections to the particle spectra due to off-equilibrium 
corrections to the ideal particle distribution function. Any bounds for the 
viscosity (at least from this paper) would have to come from the V2 spectra. 
As Lindblom [63] and earlier work by others [80] has clarified, any observable 
computed from the auxiliary fields ir^ must agree with the same observable 
generated by the physical gradients —77 (<9*u J ). When deviations are seen the 
viscous corrections can no longer be trusted. For a freezeout surface set by 
X = 4.5 the viscous corrections agree with gradients up to 2 GeV for viscosities 
as large as r]/s = 0.2 as seen in figure 3.9. It is therefore safe to use only the 
auxiliary variable when generating spectra for this particular parameter set. 
In figure 3.10 we show a summary plot of the differential elliptic flow. We 
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now show one additional curve for 77/s = 0.133 yielding (r]/p)d fl u fJ- = 0.6 
for this particular choice of freezeout surface. We believe that this choice of 
parameters is the closest physical scenario. The lower plot of figure 3.10 shows 
the measured elliptic flow as measured by the STAR collaboration [83]. We 
do not intend to make a comparison, but simply would like to keep the data 
in mind. Nevertheless since this simulation was performed with a massless 
gas which has the largest elliptic flow, it seems difficult to imagine that the 
Tj/s > 0.35 will ever fit the data even if the initial conditions are modified 
along the lines of Ref . [78] . 

Before a realistic comparison with data can be made the QGP/hadronic 
phase transition must be taken into account. In the vicinity of the phase 
transition it is possible that the shear viscosity may become very large. Also, 
a more realistic model for the hadronic gas would be the hard sphere model 
where 77 ~ This would adjust at what point the simulation freezes out and 
would therefore effect spectrum. There is most likely a finite bulk viscosity 
due to the fluctuations of the QGP and hadron concentrations in the mixed 
phase or from chemical off-equilibrium in the hadronic phase [68]. A final issue 
that should be taken into consideration is that particles of different mass could 
possibly freezeout on different surfaces. 
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Figure 3.10: (Color online) Top: Summary plot showing t>2 for massless par- 
ticles for simulations using ideal hydro and f]/s = 0.05, 0.2. Bottom: Charged 
hadron v-i data using the standard reaction plane method as measured in Au- 
Au collisions at yjs = 200 GeV for a centrality selection of 10% to 20% [83]. 
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Chapter 4 

Dilepton production from a 
viscous QGP 

4.1 Introduction 

There is a general consensus that the early matter produced at RHIC behaves 
as a near perfect fluid [84]. This conclusion was born out of the success of 
ideal hydrodynamic descriptions [35, 85] of both hadron transverse momentum 
spectra and elliptic flow measurements up to 1.5-2 GeV/c. Although it is too 
early to draw any definitive conclusions most likely the deviations from ideal 
hydrodynamic behavior can be ascribed to dissipative effects. This has already 
been suggested in some of the recent works on dissipative hydrodynamics [9, 
86, 87, 88, 89, 90, 91, 92]. 

In addition to hadronic observables which interact strongly and therefore 
depend only on the final state of the medium, electromagnetic probes are emit- 
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ted throughout the entire space-time evolution reaching the detector without 
any final state interactions. In terms of computing observables there is a big 
difference, since the resulting transverse momentum and elliptic flow spectra 
depend only on the final freezeout hypersurface, whereas the resulting dilep- 
ton yields depend on the full space-time volume. A consistent description of 
heavy-ion phenomenology should use the same space-time evolution for both 
hadronic spectra and dilepton observables. 

In this work we calculate the first viscous correction to dilepton emission 
from quark anti-quark annihilation in a dissipative medium. The kinematic 
region when a thermal description is reliable is found by requiring that the 
viscous corrections are small. When the viscous corrections become large a 
kinetic description is really required. The viscous rates are then integrated 
over the space-time history of a hydrodynamic simulation of RHIC collisions. 
We show how shear viscosity modifies the transverse momentum and invariant 
mass spectrum. We find that the inverse slope of the transverse mass spectrum 
is sensitive to both the thermalization time as well as the shear viscosity and 
can therefore be used in order to learn about the early stages of heavy-ion 
collisions. Finally a comparison is made with dileptons produced from a free- 
streaming quark-gluon plasma. 
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4.2 Dilepton Rates 



The rate of dilepton emission from a quark-gluon plasma due to qq annihilation 
is given in the Born approximation as 

Sf = / Si^iW.^/^.^^h -*■-*>)• (4.i) 

where q = (g , q) is the virtual photon's four momentum and M 2 = (E\ + 
E 2 ) 2 — (ki + k 2 ) 2 is the photon's invariant mass. Throughout this work we 
consider massless quarks; therefore = \J^\ 2 + m q ~ | ki,2 1 - The function 
f(E,T) is the quark or anti-quark momentum distribution function, which in 
thermal equilibrium is given by f(E,T) = 1/(1 + e E ^ T ). 

In the above expression V12 is the relative velocity of a quark anti-quark 
pair and a(M 2 ) is the qq cross section. Both expressions are well known from 
the literature [93] and are given by Vu = 2 e*e anc ^ ^(M 2 ) = —Jj^jt— 1 - The 
integral over the quarks' momentum can be done analytically with the result 




where n± = i/( e ('?o±|q|)/2T + ^ and f b ^T) = l/(e q °/ T - 1). 

4.3 Viscous Correction to the Dilepton Rates 

In order to account for dissipative effects in the dilepton emission rate we in- 
clude the first viscous correction to the quark and anti-quark's distribution 
function in eq. 4.1. This approach neglects any space-time inhomogeneities 
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and assumes that the distribution functions relax to their dissipative forms 
much quicker than the medium evolves. Ideally, one could solve the Baym- 
Kadanoff equations which would take non-equilibrium effects into account. We 
note that the leading order born qq rates do not contain pinch singularities 
which suggests that at least to leading order one may neglect space-time in- 
homogeneities [94]. This approximation allows us to calculate the dilepton 
emission rates locally in a space-time volume d A x. 

As shown in [51, 95, 96] viscosity modifies the ideal distribution function. 
The resulting correction for fermions is 



/(*) - /(*) + 2 {e + \)T* m[l " f(k)\k a k^ , (4.3) 



where 7r Q/ g = r)(V ' a up) = r]{S7 a up + V pu a — |A Q , j gV p M p ) and t] is the shear 
viscosity not be confused with the space-time rapidity r] s . The coefficient 
C\ can be computed analytically for a massless fermion gas and is given by 
C\ = 147r 4 /1350C(5) ~ 0.97. Substituting the above result into the born 
annihilation rate (eq. 4.1) and keeping terms up to first order in r]/s (quadratic 
in momentum) one obtains: 



dN 4:N c a 2 e 2 q f d 3 kid 3 k 2 



d A q 3(2tt) 5 J E X E. 



5 (q — ki- k 2 ) x 



>2 



f(h)f(k 2 ) + ( —^— f{ kl )[l - f(h)]f(k 2 )k^ + k 2 ^j vr Q/3 



AN c a 2 ei f d 3 kid 3 k 2 



5 (q -h- k 2 ) x 



3(2tt) 5 J E X E 2 

(4.4) 
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In simplifying the above result we have used the fact that the permutation 
of k\ <-»• k 2 gives the same contribution after integration. We write the final 
result as the sum of the ideal and viscous correction 



with 

2T ,n + . 
1 + — ln(-t) , 
|q| u- 

If = ~^yT J z^f(E 1 )[l-f(E 1 )]f(E 2 )k<?k P 1 5(E 1 + E 2 -q ). 

(4.6) 

Since if is a second rank tensor depending only on u a and q a it can be 
decomposed as 

if = a g al3 + ai«V + a 2 q a q p + a 3 («V + uV) • (4-7) 

The final result will contain the term lfit a p. By making use of the iden- 
tities u a 7r a /3 = g al3 7r a /3 = only the term with coefficient a 2 will be non- 
vanishing. a 2 is found by using the identity a 2 = P a plf where the projection 
operator in the local rest frame of the medium is 

P a p = ttm [( 3 <?o - |q| 2 )^/3 - 6qoU a qf3 + 3q a q p + \q\ 2 g a/ 3] . (4.8) 

We now quote the final result for the first viscous correction to the born 
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dilepton annihilation rates: 



dN 




v I n 

|q] n_ 



Ci 



12vr 4 



2(e + p)T 2 



(4.9) 



where we have defined 



h(qo, |q|) 



X 




and E± = |(go ± |q|)- For large invariant masses (M/T ^> 1) one can replace 
the Fermi distribution with the classical Maxwell-Boltzmann distribution. In 
the viscous correction to the distribution function this amounts to substituting 
ff{l — ff) — > Jmb- In this limit an analytic expression can be found for the 
viscous correction to the dilepton rates. In the limit that (u ■ q)/T ^> 1 the 
resulting expression is given as 



where as before C\ ~ 0.97. We find that the above result holds at the accuracy 
of a few percent for M > 3 GeV at T = 400 MeV. 

A feature of the viscous correction is that it does not modify the invariant 
mass spectrum. This is seen by looking at either of the above forms of the 
viscous correction (eq. 4.11 or 4.9). In going from d 4 q to dM 2 the resulting 
integral will be a second rank tensor depending on u a only. The most general 
form the result can take is a linear combination of terms proportional to g a/3 



dN _ N c a 2 e 2 q 
A ~ 12vr 4 



l e -1o/T 1 + 



3(e + p)T 2 



a a 



(4.11) 
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and u a u^ ', which both vanish when contracted with ir a p. 

Before performing the full space-time evolution we illustrate the effect of 
the viscous correction using a simple model for the gradients. We consider a 
ID Bjorken expansion without transverse flow. The viscous component of the 
stress-energy tensor can be easily computed and is given as 



q a q P (V a u p ) = A g 2 _ ±2 sinh 2 {y _ Vs) (4 12) 

or or 



By substituting the above result into eq. 4.11 an analytic expression can 
be found for the dilepton yields in the limit that M/T ^> 1. After performing 
the integration over rj s the result is 



dN N c a 2 e 2 q ( 1C X 



^<4+i^)[(i) 2 - 2 (^)ii) 



dM 2 dq\dy 12tt 3 

(4.13) 

where K v (x) is the modified Bessel function evaluated at x = m±/T. Fig. 4.1 
shows the kinematic regions where the viscous correction is small 
(i.e. dN vis / dN Meal < 0.8) and therefore dictates when using a thermal de- 
scription of dilepton production is suitable. The criterion that dictates when 
hydrodynamics is applicable can be written as 



2) x ± « 1 (4.14) 



and can therefore be separated into a condition on the medium, t]/s, and a 
condition on the experimental setup, 1 / (rT). Throughout this work we always 
set rj/s = 0.2. The region surrounded by the solid line is for l/(rT) w 2.2 
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Figure 4.1: Kinematic regions where the viscous correction is less then order 
one. More precisely, the boundary is set by the condition \5f/fo\ < 0.8. 

corresponding to a temperature of 450 MeV at r = 0.2 fm/c. The region 
surrounded by the dotted line is for 1/(tT) ~ 0.65 corresponding to a tem- 
perature of 300 MeV at r = 1 fm/c. At earlier times the viscous correction is 
larger and the allowed region is smaller. 

The results shown in fig. 4.1 should only be taken qualitatively. Transverse 
flow alleviates the situation, opening up the boundaries shown above. Also, 
the result presented was in the limit M/T » 1 where an analytical result 
was obtained. Fig. 4.1 is still useful, since it still serves as a qualitative pic- 
ture where the viscous corrections become large even after including flow and 
relaxing the classical limit. Outside of the kinematic boundaries a thermal 
description may no longer be reliable. 
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Throughout the remainder of this work we now resort to eq. 4.9, which is 
evaluated numerically, in order to compute the dilepton yields accurately at 
all masses. Fig. 4.2 shows the dilepton spectrum generated for a temperature 
T=0.4 GeV at proper time r = 1 fm/c and using a viscosity to entropy ratio 
of 7]/s = 0.2. 

Figures of the invariant mass spectrum are not shown because, as discussed 
earlier, the spectrum is unmodified when including the viscous correction. 
Looking at the transverse momentum spectrum a hardening of dileptons is 
seen that is reminiscent of the single particle spectrum in [51]. The magnitude 
of the viscous correction is dictated by both 7]/s as well as the proper time. 
At earlier times the shear between the longitudinal and transverse directions 
is larger resulting in bigger corrections at smaller proper times due to the 1/r 
factor in eq. 4.12. 

4.4 Evolution Model 

In order to model the space time evolution of the collision region we use the 
results of [9], which is summarized in this section. The hydrodynamic model 
is a 1+1 dimensional boost invariant expansion with initial conditions tuned 
in order to simulate Au-Au collisions at RHIC energies (-^/i = 200 GeV). 
Dissipative corrections to the ideal hydrodynamic expansion is treated using 
a second order relaxation scheme first proposed by Ottinger and Grmela [60]. 

The hydrodynamical model uses an ideal gas equation of state p = |e. 
The relationship between energy density and temperature is chosen in order 
to mimic the entire phase region through the mixed phase of an ideal Nf = 3 
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Figure 4.2: (Color online) Dilepton transverse momentum spectra for T = 0.4 
GeV and f]/s = 0.2 at r = 1 fm/c for a boost invariant expansion with no 
transverse flow. 
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QGP equation of state with a 1 st order phase transition. This happens to 
correspond to an ideal gas having an effective DoF pa 16. At an initial time 
the entropy is distributed in the transverse plane according to the distribution 
of participants for a Au-Au collision. Then one parameter, C s , is adjusted to 
set the initial temperature and total particle yield. The value C s = 15 closely 
corresponds to the results of full hydrodynamic simulations and corresponds 
to an initial temperature of To = 420 MeV for tq = 1 fm/c. 

Detailed plots of the hydrodynamic solution with and without viscosity is 
shown in [9]. For modest values of the shear viscosity (rj/s ^ 0.3) dissipative 
effects did not integrate to give large changes to the ideal hydrodynamic so- 
lution. The net effect of a finite viscosity was twofold. First, the longitudinal 
pressure is initially reduced causing a slower decrease of energy density per 
unit rapidity at early times. The reduction of longitudinal pressure is accom- 
panied by a larger transverse pressure which drives larger transverse velocities. 
The larger velocities then cause the energy density to deplete faster at later 
times. 

Even though the changes to the ideal hydrodynamic result is small a full 
viscous simulation is still needed in order to have access to the velocity gradi- 
ents which enter into the dissipative corrections of the quark and anti-quark 
distribution functions. 

The hydrodynamic model is started at To = 0.2 fm/c in order to account 
for some of the pre-equilibrium production of dileptons which will contribute 
at larger masses. Dileptons are produced as long as the temperature of the 
medium is greater than a critical temperature taken as T c = 0.170 GeV. We do 
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not look at dileptons produced during a mixed phase or hadronic phase in this 
work. At any space-time point the hydrodynamic model provides the three 
terms of the stress tensor (7T rr ,r 2 iT^ and r 2 7r w ). Actually for an azimuthaly 
symmetric hydrodynamic simulation only two terms of the viscous stress tensor 
are independent since the third term could have been found by making use 
of the tracelessness of the shear tensor as well as its orthogonality to [97]. 
The equation for q a q /3 ir a p used in eq. 4.9 is given as 

q a qPit a p = qlcos\9)n rr + qlsm\9)r\^ + mlsmh 2 {r ]s )T 2 ir rm 

+m 2 L cosh 2 (?7 s )w 2 7r rT — 2m± cosh(r] s )q± cos(9)vir rr , (4.15) 

where 9 = <p q — <p v is the relative angle between the virtual photon's momentum 
(q±) and the fluid cell's radial velocity. 

Fig. 4.3 shows the resulting transverse momentum spectrum after the full 
space-time integration at two invariant mass points: M = 0.525 GeV (top) 
and M = 2.625 GeV (bottom). First the red curve shows spectra generated 
from an ideal hydrodynamic simulation (r]/s = 0). Next the green curve shows 
the spectra generated from a viscous simulation having r]/s = 0.2 but without 
including the viscous correction to the distribution function. This curve there- 
fore shows the effect that viscosity has on modifying the ideal hydrodynamic 
equation of motions. We find that a finite viscosity leads to a slight increase 
in the overall yield. This is due to the fact that a finite viscosity causes the 
energy density to deplete more slowly at early times. This effect therefore 
brings about an effective increase in the lifetime of the simulation above the 
critical temperature. We find ~ 30% increase in the low mass region and 
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~ 50% increase in the higher mass region. 

Finally, the blue curves in fig. 4.3 show the viscous result including the 
viscous correction to the distribution function. We find that the magnitude of 
the viscous correction increases with the invariant mass. This was similarly 
observed in fig. 4.1 where the range in q± having viscous corrections of order 
less than one (as shown by shaded regions) decreased in size with increasing 
mass. The simulation results are discussed in more detail in the next section. 



4.5 Discussion 

In order to further understand the viscous corrections, the effective temper- 
ature (T e ff) of the dilepton spectrum from the full space-time integration is 
constructed as a function of invariant mass. The effective temperature is found 
by fitting the transverse mass spectrum at a given mass to the following ex- 
pression, 

oc e - m ^ff . (4.16) 



dM 2 m±dm±dy 

In this work the fit is done in the transverse momentum region of 0.5 < 
g_i_(GeV) < 2.0. As expected, we find that the transverse mass spectra does 
not exactly fit the above form. Actually, other ranges in q± could have been 
chosen where the fit works better. However, the results are qualitatively the 
same and therefore a different choice in range will not change the discussion 
that follows. If a quantitative comparison were to be made with data, it would 
be more appropriate to compare to the actual q± spectra instead. Regardless, 
T e ff still serves as a useful quantity since it probes the average temperature of 
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Figure 4.3: (Color online) Dilepton transverse momentum spectra after the full 
space-time integration of a boost invariant expansion with arbitrary transverse 
expansion and azimuthal symmetry. 
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Figure 4.4: (Color online) Effective temperature as a function of invariant 
mass. 

the medium as well as the radial flow profile and viscous correction. Looking at 
fig. 4.3 we expect the viscous correction to increase the effective temperature, 
with larger corrections at higher masses. 

In fig. 4.4 the effective temperature is shown as a function of invariant 
mass. The solid blue curve labeled ideal shows T e ff for an ideal (rj/s = 0) 
hydrodynamic expansion started at proper time tq = 0.2 fm/c. The shape 
of the curve is dictated by the underlying radial flow as well as the average 
temperature of the emission region. At higher masses the dominate source of 
dileptons is from the higher temperature regions, which occur at earlier proper 
times. This explains the slight rise in T e ff with mass. We now look at the solid 
red curve, which is generated from a hydrodynamic evolution having rj/s — 
0.2. In this case we do not include the viscous correction to the distribution 
function and the resulting modifications to the effective temperature are due 
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to changes in the hydro dynamic evolution. As discussed earlier, modifications 
from viscosity to the hydrodynamic solution are small and we therefore don't 
expect to see large deviations from the ideal case. This is indeed the case. 

We now focus the discussion on the role of the viscous correction to the 
distribution function. The green curve in fig. 4.4 shows the effective temper- 
ature of dileptons coming from a viscous medium having rj/s = 0.2 from a 
simulation started at a proper time of tq = 0.2 fm/c. The result is that the 
effective temperature increases greatly as a function of invariant mass. From 
the magnitude of the correction, the upper bound of the domain of hydrody- 
namics is found to be at most M max m 2.0 GeV for this parameter set. There 
are two reasons why the viscous correction increases with mass. First, there is 
the explicit mass dependence in the viscous correction itself. This is easiest to 
see by looking at the approximate form, eq. 4.13. The second reason is because 
the high mass contribution is mainly produced in the early, high temperature 
stages of the evolution. Looking at eq. 4.12 the viscous correction grows like 
1/r at early times. In order to see the effect of the early emission a final 
simulation is done (dashed-black curve) where the hydrodynamic evolution is 
started at To = 1 fm/c. In this case the viscous corrections are more modest 
and M max ^4.5 GeV. 

It is therefore a combination of both the thermalization time as well as the 
magnitude of rj/s that dictates when a hydrodynamic description is reliable. 
Since the effective temperature rises so quickly with mass, as long as there is 
non- vanishing viscosity, there will always be a mass region outside of the region 
of a hydrodynamic description. From eq. 4.13 one can extract an approximate 
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condition for the mass. Since most of the particle yield is at low q± we should 
guarantee that the viscous correction is small at q± = 0. Furthermore at high 
and intermediate masses the ratio of Bessel functions is approximately one. 
Then the condition that the viscous correction must be less than of order one 
can be expressed as 

M max » ^ . (4.17) 
rj/s 

When the viscous corrections to the spectra become large a kinetic ap- 
proach is required. One can ask whether the viscous correction at early times 
mock up the effects of off-equilibrium production that would be taken into 
account by a full kinetic theory. In order to test this hypothesis T e ff spectra 
is calculated from a free streaming non-interacting gas of quarks [98, 99]. We 
should point out that our treatment is very similar to the recent work of [100]. 
In this model the initial parton distribution is taken as thermal with the tem- 
perature chosen in order to reproduce the thermal dilepton number given by 
the hydrodynamic simulation. Starting with the thermal initial condition at 
r = 0.2 fm/c the total dilepton yield is found by integrating the free streaming 
result 4.27 up to a final time of r = 1.0 fm/c. The details of this calculation 
is given in the appendix. We now discuss the results. 

We consider two scenarios. The first is running the hydrodynamic simu- 
lation starting at r = 0.2 fm/c until T c . The second scenario runs the free 
streaming model from 0.2 < r(fm/c) < 1.0. Then at 1.0 fm/c the hydrody- 
namic evolution is started and ran until T c . We should stress that the second 
model is not very realistic since the free streaming model is not asymptotic 
with the hydrodynamic evolution at r = 1 fm/c. A future work might use 
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a more realistic model for the evolution of the distribution function then the 
proposed free streaming case. For example, one could start with an initially 
anisotropic distribution which evolves to its thermal form from multi-quark 
scattering [101], at which point a hydrodynamic evolution is started. How- 
ever, we expect the true result to lie between the two scenarios used in this 
work. 

Figure 4.5 shows both the invariant mass spectrum (top) and effective 
temperature (bottom) for the two scenarios outlined above. We note that the 
invariant mass spectrum is generated by integrating over all q±. The curves 
to compare are the hydrodynamic simulation started at 7~o = 0.2 (labeled Hy. 
To = 0.2 fm/c) and the sum of the hydrodynamic simulation started at r = 1 
and the free-streaming (labeled FS+Hy.). We first note that the invariant 
mass spectrum is qualitatively the same for the two scenarios and it would 
not be possible to discern between the two scenarios from experimental data. 
Qualitative differences do appear in the T e ff spectrum. First we find that the 
free-streaming with hydro solution mimics the early time hydro only solution 
at low masses. However, at high masses the two result diverge when the viscous 
correction can no longer be trusted. While the early hydrodynamic solution 
increases greatly with mass the free streaming solution flattens off at higher 
mass. We therefore argue that through a detailed analysis of q± spectra one 
could hopefully extract information on the thermalization time, viscosity to 
entropy ratio and thermalization mechanism in heavy-ion collisions. 

Some comments are in order on how viscosity may modify the dilepton 
emission from the hadronic phase, which comprises of a larger part of the 
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Figure 4.5: (Color online) Top: Dilepton invariant mass spectra. Bottom: 
Effective temperature as a function of invariant mass. 
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overall yield in the mass region considered here. A more detailed analysis will 
be presented in a future work. As the medium evolves the viscous pressure 
decreases rapidly so one may expect dissipative effects to be smaller. However 
the viscosity becomes larger as the temperature decreases. For a pion gas in 
the chiral limit f]/s = jjr- (^) 4 and it therefore becomes a dynamical question 
on how large the viscous corrections become in the hadronic phase. 

4.6 Conclusions 

In conclusion, we have calculated the first viscous correction to dilepton pro- 
duction from leading order qq annihilation. The rates are then integrated 
over the space-time history of a viscous hydrodynamic simulation of RHIC 
collisions. We argue that a thermal description is only reliable for invariant 
masses less than m (2tqTq) / (rj / s) and above this a kinetic description is re- 
quired. For rj/s — 0.2 and r = 1 fm/c this corresponds to M < 4.5 GeV. 
We have shown that viscosity does not change the invariant mass distribution 
but strongly modifies the transverse momentum distribution and can therefore 
be used to extract information on both the viscosity to entropy ratio as well 
as the thermalization time. Finally, we have also made comparisons with an 
initially free streaming QGP. Qualitative differences in transverse momentum 
are seen, which could again possibly be used to learn about the thermalization 
mechanism. 



92 



Appendix: Dilepton yields from a free stream- 
ing QGP 



In this appendix the dilepton yields are derived for a free streaming boost in- 
variant expansion. The starting point is the collision-less Boltzmann equation 

jfd ll f( P ,x) = 0, (4.18) 

where f(p,x) will be considered as the phase-space distribution for the quark 
and anti-quark. Under the assumption of boost invariance as well as homo- 
geneity in the transverse plane the Boltzmann equation can be written as 



Orf-^^dJ-O, (4.19) 



where x = V ~ Is- The initial condition is such that the quark distribution 
is isotropic and starts from local thermal equilibrium, f(p,r = r ) = eP0 /T +1 • 
One can write Po = u-p = pj_ cosh(x) by using the assumption of boost invari- 
ance and homogeneity in the transverse plane. We note that even for quarks 
out of equilibrium it still might be useful to use the equilibrium form of the 
distribution function where T is instead considered as an effective temperature 
describing the initial state. The solution of eq. 4.19 at any time r is 

f(P, x) = 1 . (4.20) 

e ^V 1+sinh2(x) te) +i 

With the explicit form of the distribution function available one can cal- 
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culate the dilepton rates using the same kinetic theory expression used before 
(see eq. 4.1) 



dN 
~¥q 



Making use of the expressions for the relative velocity and cross section as 
quoted earlier the above equation can be expressed in the following form 



B I d'x I d%d%5(pl)5(pl)5^( Pl +p 2 -q)f(p l ,x)f(p 2 ,x), (4.22) 



d A q 



where B 



First let us quote some well known identities: 

d A x = TdTdrj s d 2 x± = irR 2 TdTdr] s 

d 4 p = -dp 2 dy p pj_dp±d(p p 

d 4 q = MdMq ± dq ± dyd(j) 

6W{P-q) = 46(P 2 -M 2 )6(y p -y)6(P L -q ± )5( ( p P -^ 

where P M = P1+P2 and y and are rapidity and angle in the transverse plane. 
We place a subscript P on quantities to indicate they are derived from P M . 
The free streaming dilepton rate can now be expressed as 



B I d 4 x I dy 1 p 1)± dp 1)± d(j) 1 dy 2 p2,±dp2 ) ±d(j)2f(pi,x)f(p2,x) x 



d A q 



5(P 2 - M 2 )5(Pi - q 2 ± )5(y P - y)5(<j> p - 0) (4.23) 
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Since the distribution function is boost invariant the integral over rj s is 
trivial due to the delta function, 5(yp — y). By defining y± = yi ± y 2 and 
<j)± — 0i ± (f> 2 the delta functions can be rewritten as 



5(P 2 - M 2 ) = - 5(cosh y_ - cos ' 



2pi,_LP2,_L 2Pl,±P2, 

2 i J2 Jl 



X(r>2 2\ 1 r / / I Pl,_L+P2,± <?±n 

<H^±-<?J = ^ (5(COS0_ + ) 

After rewriting the integration variables as dy\dy 2 = dy+dy ^ and dcpidfo 
dc/> + d<i>~ ^ e j n ^ e g ra | over 0_ anc | ^_ can ci orie explicitly yielding 



dN - 2d2d 



MdMdyq ± dq ± 

Pi,±dpi,±p 2 ,±dp 2 , 



J rdr J dy+ (2p h± p 2: ±) 2 | sin 0_ | sinhy„ 
x f(pi,x)f(p 2 ,x) (4.24) 



where 

i • , i V ((Pi,± + P2,±) 2 - oY) (Q± ~ (Pi,-L - P2,±) 2 ) 
sm0_ = — 

2Pl,±P2,X 

ginh y/(M 2 + q{ - ( Pl , ± +p 2 ,±mMi + q{- (p h± - p^f) 

%Pl,XP2,X 

(4.25) 



The delta function in the above equation enforces the following constraints 



Pi 



pI 



2Pl,±P2,± 

M 2 + q\- p 2 1± - pj L 

2>Pl,±P2,± 



^ 1 
^ 1 



(4.26) 
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Let us make a further shift of variables, p± = pi t ± ± p2,±- The constraints 
(4.26) then take particularly simple form. The final expression is 



dN 



dM 2 dydq] 



,N c a 2 e 2 q 

487T 4 



+00 



rdr / dy^ 



g± 



dp + I dp. 
-g± 
1 



V / (M 2 + qi ~ f+){M 2 + qi _ pl) _ qi){qi _ pi) 

x (pl-pl)f(Pur)f(p 2 ,r) (4.27) 



where 



/(Pi.r) 



f(p2,r) 



p + — p^ 1 1 r \ I y+ — U- 
1 + exp| ^^\/ 1+ Uj sinh l^ 



sinh 1 



2y/(M 2 + qj- gj (M 2 + q\ £g 



(4.2? 
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Chapter 5 

Heavy Ion Phenomenology 



5.1 Di-muons at the CERN SPS collider 
5.1.1 Introduction 

It is expected that above a critical temperature, T c 170 MeV, QCD under- 
goes a chiral phase transition where the relevant degrees of freedom change 
from mesons and baryons to a phase of strongly coupled quarks and gluons, 
the strongly coupled quark-gluon plasma (sQGP). This new phase of matter is 
being searched for in a number of ultra-relativistic heavy ion facilities. There 
are a number of current observations in favor of the sQGP, ranging from hydro- 
dynamical flow (soft probes) to jet quenching (hard probes). However, most 
of these observations are blurred by the fact that the space-time evolution of 
the sQGP is short and its conversion to hadronic matter is involved. Since 
the latter dominates the final stage of the evolution, it is producing competing 
signals that interfere with those from the sQGP. In this respect, dilepton and 
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photon emissions are interesting probes of the collision region as neither inter- 
act strongly with the medium produced in these collisions, thus they probe the 
early stages of the collision. This is in contrast to hadronic observables which 
thermalize along with the collision region thus providing information only on 
the late (or freeze-out) stage of the collision. 

Making quantitative predictions of the production rates of dileptons and 
photons is difficult for a number or reasons. Since the temperature produced 
in typical heavy-ion collisions is in the range of 200-300 MeV which is about 
the QCD scale factor, A, the differential cross sections can not be computed in 
a weak-coupling expansion. Another uncertainty is detailed knowledge of the 
evolution of both hadronic matter and quark gluon phase produced in heavy- 
ion collisions. In addition there is also a background of dileptons from other 
processes not occurring in the quark-gluon plasma such as hadronic decays. 

In the past there have been a number of experiments probing photons and 
dileptons created in hadronic collisions. One of the most recent experiments 
was the CERES (NA45) taking place as the CERN SPS collider which looked 
for dielectrons. It was found that the dielectron production exceeded the 
theoretical expectations for conventional processes in both hadronic and QGP 
matter [102], especially in the mass region 0.3 < M(GeV) < 0.6 [103]. A 
number of theoretical analyses were put forward to explain this excess based 
on effective Lagrangians with medium modification [104, 105] and dropping 
vector meson masses [106]. Model independent emission rates constrained by 
the strictures of broken chiral symmetry and data were unable to account for 
the excess rate reported by NA45 [6, 7, 8]. However, the large statistical and 
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systematic errors reported by NA45 in exactly the excess region, did not allow 
for a definitive conclusion as to the theoretical nature of the emissivities. 

In this letter we revisit these issues in light of the recently reported dimuon 
data from the NA60 collaboration using In-In collisions at 158 Gev/Nuc [107]. 
These data have far better statistics, which gives much better constraints 
on any medium modification to the vector mesons [108, 109]. We use the 
model independent analysis in [6, 7, 8] to analyze these data, whereby the 
emissivities are constrained by broken chiral symmetry in a dilute hadronic 
medium, and by non-perturbative QCD in the sQGP. The collision expansion 
and composition are extracted from an underlying hydro dynamical evolution 
set to reproduce the CERN SPS conditions. 



5.1.2 Dilepton Emission Rates 

The rate of dilepton emission per unit four volume for particles in thermal 
equilibrium at a temperature T is related to the thermal expectation value of 
the electromagnetic current-current correlation function [21, 22]. For massless 
leptons with momenta pi and P2, the rate per unit invariant momentum q = 
Pi + P2 is given by: 

7 q = S$TT^ ImWFW (5 - 1) 

where a = e 2 /4ir, T is the temperature and 

W F (q) =tjd 4 x e iqx Tr [e- {n '^ N - n)/T T*3^(x) J M (0)] (5.2) 
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where eJ M is the hadronic part of the electromagnetic current, H is the hadronic 
Hamiltonian, \xb is the baryon chemical potential, N is the baryon number op- 
erator, and Q is the Gibbs energy. The trace is over a complete set of hadron 
states. 

In order to take into account leptons with mass mi the right-hand side of 
Eq. 5.1 is multiplied by 

2m? , , 4m? , , /9 . . 

(1 + ^)(1-^) 1/2 (5.3) 

To compare the theoretical dilepton production rates with those observed in 
heavy ion collisions, the rates must be integrated over the space-time history 
of the collision region and then finally integrated over the dilepton pair's trans- 
verse momentum and rapidity in order to compare with the yields measured 
by the NA60 collaboration. The final expression for the rates is given as: 



dN 
dM 



2,Mj d yJ dq± . qx ,A CC (M, q± ,y) 

d 4 qd A x 



fTf pco pr ma x r2n d^R 

x / rdr / drj rdr / dO (M, \q\,T, fi B , x) 



o 



(5.4) 



where M = \[q* is the dilepton invariant mass, y is the dilepton pair rapidity, 
rj is the spatial rapidity, q± is the dilepton pair transverse momentum (with 
defined as the angle between q± and the fluid element's velocity), x is the 
hadron fraction, r is radial coordinate (with r max set by the freeze-out temper- 
ature), and Acc(M, q±,y) is the experimental acceptance taking into account 
that the CERES detector covers a limited rapidity in the interval y = 2.9 — 4.5 
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in the lab frame. 

The integration over rj, r, 9 and r was done over the full hydrodynamic 
simulation of the collision region as described below. \q\ can be found by 
considering the two invariants q^q^ = M 2 and u^q^ constructed from the 
dilepton momentum and fluid 4-velocity which can be expressed as: 

g M = (M ± cosh(y), q ± , M ± sinh(?/)) 

u 11 = (7± cosh(?7), 7^wj_, 7j_ sinh(ry)) (5.5) 

giving 

\q\ = [-M 2 + (7 ± M ± coshfo) - u ± q ± cos(#)) 2 ] 1/2 (5.6) 

where u± = 7j_w±, 7_l = , 1 „ , and is the transverse fluid velocity which 
is taken from the hydrodynamic simulation. 

The acceptance function has a complicated dependence on M, q±_ and y, but 
since our rates are y-independent we have used an acceptance with M and q± 
dependence built to specifications provided by the NA60 collaboration [110]. 
Without detailed hadronic data available (such as transverse mass spectra 
and HBT analysis) a careful consideration of hadronic input, such as freeze- 
out temperature, cannot be made. Therefore there is a large uncertainty in 
the overall normalization of the yields, which depends strongly on T/. .. In 
addition, the assumption of boost invariance can also affect the normalization 
as the acceptance probes very specific rapidities. The approach taken here is 
to normalize our results to the excess data in the peripheral centrality windows 
which fixes the normalization in the central bins. 
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5.1.3 Spectrum above Tc 

At temperatures T>Tc lattice calculations have predicted that the relevant 
degrees of freedom consists of (strongly) interacting quarks and gluons. In 
order to compute the dilepton production rates as one would expect from a 
conventional phase of quark-gluon plasma we use the Born qq annihilation 
term [111, 112] which for massless quarks is 



\ q=u,d,s / 



\q\ n. 



(5.7) 



where Nc is the number of colors, e q the charge of the quarks, and n± = 
l/( e (qo±\<i\)/ 2T _|_ i). it should be mentioned that 5.7 reduces to the well-known 
vacuum result (ImW^ = ~ q 4 ^ c Y2 q e q) a ^ T=0 and that the finite temperature 
rate is always smaller then the T=0 rate due to Pauli blocking. 

It has also been seen in lattice simulations that near the critical temper- 
ature T c there are still substantial chromoelectric and chromomagnetic con- 
densates present leading to additional non-perturbative effects. It was shown 
in [113] that the enhancement to the dilepton rates in a plasma with non- 
vanishing chromoelectric and chromomagnetic condensates can be given by 



ImW R = 

4:71 



\ q=u,d,s / 



q 2 ( ^ A 2\_ 1/^1/^2 

7i I 6\7r / 3 \ 7i 



Att 2 , 

(n+(l — n + ) — n_(l — n_)) (51 



T\q 



where the values of the above condensates in 5.8 can only be estimated us- 
ing non-perturbative calculations such as lattice QCD. The net result is a 
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Figure 5.1: Integrated dimuon rates from the plasma phase for T=150 MeV. 
The thick solid line shows the perturbative qq annihilation rates while the 
this solid lines show the results for non- vanishing A 2 , B 2 and E 2 condensates 
and for only a non- vanishing A 2 condensate. For comparison the integrated 
hadronic rate at T=150 MeV and /is = 225 MeV is also shown, which will be 
discussed in the next section. 

substantial enhancement (as seen in Fig. 5.1 for the case of dimuons) of the 
dilepton production rates below an invariant mass of ~ 0.4GeV. For the re- 
mainder of the paper we refer to the perturbative plasma rates as those given 
by 5.7 and the non-perturbative plasma rates as those given the sum of equa- 
tions 5.7 and 5.8 using 5 = (^E 2 ) /(200MeV) 4 = (^5 2 ) /(200Me1/) 4 = 
and (3 = (^A 2 ) /T 2 = 0.4 which is the lower of the non-perturbative curves 
in Fig. 5.1. An explanation of the choice of f3 = 0.4 can be found in [113]. 
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5.1.4 Spectrum below 

Even though there are various approaches to calculating production rates, they 
differ in the way in which the current- current correlation function in Eq. 5.1 
is approximated and evaluated. The approach taken here is to use a chiral 
reduction formalism in order to reduce the current- current correlation function 
in 5.2 into a number of vacuum correlation functions which can be constrained 
to experimental e + e~ annihilation, r-decay, two-photon fusion reaction, and 
pion radiative decay experimental data. 

For temperatures, T < and for nucleon densities, piv < 3po the trace 
in Eq. (5.2) can be expanded in pion and nucleon states. Keeping terms up to 
first order in meson and nucleon density gives [8] 

ImW F (g) = -3g 2 ImIIv,(g 2 ) + J daWf(q, k) + J dNW F N {q,p) (5.9) 

with phase space factors of 



d 3 v 1 1 
dN=^-— — r= =-= (5.10) 

(2tt) 3 2E p e ( E p~^)/T + 1 v J 



and 



d 3 k 1 1 . 

da = - — — -rj= (5.11 

(2tt) 3 2ul e^ T - 1 V ' 



with nucleon and meson energies of E p = a/ m 2 + p 2 and uj^ = \Jm 2 a + k 2 
respectively. 

The first term in 5.9 is the transverse part of the isovector correlator 
(0|T*VV|0) which can be determined experimentally from electroproduction 
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data and gives a result analogous to the resonant gas model. At low and in- 
termediate invariant mass the spectrum is dominated by the p(770 MeV) and 
p'(U50 MeV). 

The term linear in meson density (the second term in Eq. 5.9) can be related 
to experimentally measured quantities via the three flavor chiral reduction 
formulas [23]. It is shown in [6, 8] that the dominant contribution comes 
solely from the part involving two-point correlators which gives: 



Wf (g, k) = ^g 2 Imn y (g 2 ) + -^g 2 Im (V y (g 2 ) + ^(g 2 )) 
_|_ (fc + g) 2 Iinii; ((* + g )») + (g - -g) 

- + q f [l m n^ ((* + g) 2 ) + Imn^ ((* + g) 2 )] + (g - -g) 

Jk 

+ 4 ((& • qf - ™W) Imn 7 y (g 2 ) x ReA^fc + g) + (g - -g) 

J TT 

+ -L({k. qf - mW) Im (W(g 2 ) + X ReAf (fc + g) + (g - -g) 

(5.12) 

Where is the retarded meson propagator given by l/(g 2 — m\ + ie) and 
11^ is the transverse part of the iso-axial correlator (0|T*j J 4j_A|0). The spectral 
functions appearing in Eq. (5.12) can be related to both e + e~ annihilation 
as well as r-decay data as was compiled in [24]. As already shown in [6] 
the spectral function can be directly related the form factor, Fy, via the 
KSFR relation where Fy is parameterized in the common Breit-Wigner form 
where the resonance parameters and decay constants are taken from empirical 
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Figure 5.2: The total integrated dimuon rates from a hadronic gas at T=150 
MeV. The curve labeled Res. Gas shows the analogue of the resonant gas 
contribution (the first term in Eq. 5.9). The points labeled Tly and II a give 
the contributions from all of the respective spectral functions in equations 5.9 
and 5.12. The thin line labeled meson is the total rate given by a hadronic 
gas without nucleons. The thick solid line gives the total dimuon rate when 
nucleons (shown here for \Ib — 225 MeV) are taken into account. 



data. Included in the data are contributions to the spectral function from the 
p, u, 0, ai, Ki and some of their radial excitations (see Table I in [8]). 

It can be seen in Fig. 5.2 that the term linear in meson density decreases the 
rates from the resonance gas contribution for the mass region above the two 
pion threshold. However below the two pion threshold the only contribution 
to the rates come from the II ^ terms in Eq. 5.12. This is because the axial 
spectral density is integrated over all momentum in the thermal averaging 
(Eq. 5.9), which weakens the (k + q) 2 factor in Eq. 5.12 allowing the 1/q 2 term 
in Eq. 5.1 to dominate at low q 2 . 

The final term in Eq. (5.9) which is proportional to the nucleon density is 



106 



the spin-averaged forward Compton scattering amplitude of virtual photons 
off a nucleon. Experimentally, data is only available for values of q 2 < 0, so 
while the photon rate which requires q 2 = can be determined by use of the 
optical theorem the contribution to the dilepton rates must be determined by 
chiral constraints. Broken chiral symmetry dictates uniquely the form of the 
strong interaction Lagrangian (at tree level) for spin-| particles. Perturbative 
unitarity follows from an on-shell loop-expansion in j- that enforces current 
conservation and crossing symmetry. To one-loop, the ttN contribution is 
parameter free. The large contribution of the A to the Compton amplitude 
near threshold is readily taken into account by adding it as a unitarized tree 
term to the one-loop result [7, 114]. The enhancement in the dimuon rates 
due to a non-vanishing baryon density can be seen in Fig. 5.2 where the solid 
curve shows the total dimuon spectra with an enhancement as large as a factor 
of two in the invariant mass region of 2m ^ < M(GeV) < 0.6. 

5.1.5 Hydrodynamic Evolution 

As mentioned earlier, in order to compare the theoretical dilepton production 
rates with those seen in heavy-ion collisions it is necessary to integrate these 
rates over the space time evolution of the collision region. We consider a 
region localized in space-time consisting of thermal hadronic matter acting 
as a source of particles. Equilibrium of the collision region is strictly a local 
property with different temperatures and baryon densities possible in different 
space-time domains. 

A computational hydrodynamic code was already developed by one of us 
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and it has been modified to the conditions of the SPS collider for Indium on 
Indium collisions. In this paper we only briefly outline the physics behind 
this code and show the results of the In-In collisions which has not been mod- 
eled before. For technical details regarding the hydrodynamic calculations the 
reader is referred to the prior works by one of us [115]. 

Hydrodynamics 

The hydrodynamic equations for a relativistic fluid consist of the local con- 
servation of energy and momentum, which can be written in compact form 
as d u T^ LU = 0, as well as local charge conservation c^Jf = where = 
(e + p)U^U u —pg^ u is the energy-momentum tensor with e the energy density, 
p the pressure, = v) is the proper velocity of the fluid, and Jf is any 
conserved current (e.g. isospin, strangeness and baryon number in the case of 
strong interactions). 

The same space-time evolution scenario as first proposed by Bjorken [116] 
is assumed where the equation of motion can be described by the Bjorken 
proper time r = \Jt 2 — z 2 and the spatial rapidity y = |ln|^. One of the 
main results, following from the assumption of a central-plateau structure in 
the rapidity distribution is that of boost invariance, stating that the initial 
conditions and thus the subsequent evolution of the system are invariant with 
respect to a Lorentz boost. Thus a solution at any value of y can be found by 
boosting the solution at y = to a new frame moving with velocity v = tanh(y) 
in the negative z-direction. 

With the assumption of boost invariance the equations of motion are a 
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function of the transverse coordinates and the proper time r only. After 
integrating over the transverse plane of the collision region one finds that 
(dStot/drj), (dnB/drj), and the net transverse momentum per unit rapidity are 
all conserved. 

Equation of State 

In order to solve the equations of motion as given by the vanishing of the 
divergence of the energy-momentum stress tensor one must have an Equation 
of State (EoS) relating the local values of the pressure, energy density, and 
baryon density (ng). The approach taken here is to consider an EoS with a 
variable latent heat in the e/ng plane where the following two derivatives hold 
along a path of constant ub/s: 



If the speed of sound is defined everywhere along with the entropy of one 
arc in the e, % plane the above derivatives can be integrated in order to 
determine s(e,ns). From the entropy all other thermodynamic variables, such 
as T and can be found as needed. 

We consider a hydrodynamic evolution that consists of three phases, a 
hadronic phase, a QGP phase, and a mixed phase. The hadronic phase is 
taken to be made of ideal gas mixtures of the lowest SU(3) multiplets of mesons 




(5.13) 




(5.14) 
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and baryons. All intensive thermodynamic quantities including p, e, s, and hb 
can be found as a sum of that quantity's contribution from each specie in the 
gas consisting of a simple Bose or Fermi distribution. The hadronic phase 
is taken up to a temperature of Tq < 170 MeV and has a squared speed of 
sound of approximately l/5c 2 . For temperatures above Tq only the squared 
speed of sound, c 2 s is specified. For the mixed phase it is taken almost at 
zero (c 2 = 0.05c). For the QGP phase the degrees of freedom are taken to be 
massless and the speed of sound is accordingly c 2 = l/3c. The extension of 
this analysis to the sQGP is beyond the scope of this work. 

Initial Conditions 

The initial conditions of the collision consist of setting the entropy and baryon 
density proportional to the number of participating nucleons in the transverse 
plane at some initial proper time r = 1 fm/c. Since both the entropy and 
baryon number per unit rapidity are conserved the final yields of pions and 
nucleons are proportional to the number of participants. The number of par- 
ticipants were calculated by use of a Glauber model and the initial entropy and 
baryon densities were fixed by two constants C s and C ng , which respectively 
are the entropy and net baryon number produced per unit spatial rapidity per 
participant. These constants were fixed to the conditions at the CERN SPS 
collider in order to fit the total yield of charged particles and the net yield 
of protons. Table 5.1 summarizes the input parameters used in the hydro- 
dynamic calculations. In order to address the centrality of the collision the 
impact parameter was chosen in order to reproduce the number of participants 
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Parameter 


Value 


2 

mixed 


U.Uoc 


2 

C QGP 


0.33c 


T 


1 1 u ivie v 


T/.o. 


130 MeV 


To 


1.0 fm/c 




0.0238 


C s 


8.06 


Cn B 


0.191 



Table 5.1: Parameters used in the hydrodynamic simulation of In- In collisions, 
as reported in [117]. 

5.1.6 Results for In-In Collisions at CERN SPS 

The hydrodynamic result for In-In Semi-Central collisions is shown in Fig. 5.3. 
The two thick lines labeled cq and e# represent contours of constant energy 
density showing the transition from the plasma phase to the mixed plasma 
and hadronic phase and the transition from the mixed to the purely hadronic 
phase respectively. It can be seen that the QGP phase takes up a much 
smaller space-time volume then the hadronic phase, however the rates still 
appear in the spectrum as the high temperatures in this region enhance the 
rates by an order of magnitude. The effect of nucleons depends on the baryon 
chemical potential throughout the evolution. This is plotted as a function of 
temperature for the pure hadronic phase in Fig. 5.4. 

5.1.7 Results and Discussion 

Our final dimuon yields for the four centrality windows is shown in Fig. 5.1.7 
where it is compared to the excess data measured by the NA60 collaboration. 
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Figure 5.3: The hydro dynamic solution for semi-central In-In collisions at 
the SPS collider. The thin lines show contours of constant transverse fluid 
rapidity (v± = tanh(y^)) with values of 0.1,. 02,.., 0.5. The dashed lines show 
contours of constant temperature with values of (working radially outward) 
T=200 MeV, T=150 MeV and T=135 Mev. The e Q = 1.70 GeV/fm 3 and 
en = 0.50 GeV/fm 3 contours represent the phase changes from QGP to mixed 
and from mixed to hadronic matter respectively. 
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Figure 5.4: Dependence of the baryon chemical potential, on the temper- 
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In all four figures we show the total dimuon yield, which includes contributions 
from the hadronic phase, either the perturbative or non-perturbative plasma 
phase, as well as the DD contribution as provided by the NA60 collaboration. 
For all cases we also show separately the perturbative and non-perturbative 
plasma contributions to the overall yield. It can be seen to be almost negligible 
in the peripheral data. For the central data, where there is a larger plasma 
contribution we also show curves showing separately the hadronic contribution. 

Even though it can be seen that the theoretical rates are able to describe 
most of the features of the spectrum, a number of things should be noted 
before a direct comparison is made. The rates below M=0.4 GeV should not 
be taken literally since they are obtained by saturating the total measured 
yield in that region by r\ Dalitz decays only, thereby lowering the excess close 
to the dimuon threshold. Actually, by reducing the 77 yield by 10% the data has 
much better agreement with the theory for M< 0.4 GeV. The charm decay 
data was analyzed and provided by NA60, and since the contribution from 
charm decays is not subtracted from the excess data it must added to our 
rates for comparison with experiment. The excess spectra which is shown in 
the figure is created from subtracting the cocktail (omitting the p(770)) from 
the total observed data. This would erase any uj or <p peak at the vacuum 
positions. Since our hadronic rates don't modify either the position or width 
of the uj or it can be very difficult to distinguish any residual uj or 0's from 
the cocktails'. 

It can be seen right away that the dimuon yields are reproduced in the 
peripheral centrality windows. This is expected as the matter is dilute and 
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any medium modification to the spectral densities will be accounted for in the 
virial expansion (Eq. 5.9). In the central bins it can be seen that the shape of 
the spectrum changes as one goes to more central collisions. Even though the 
general shape of the spectrum is reproduced by our rates, our rates slightly 
over-predict the yield at the p peak by about ~ 50% for semi-central and by 
~ 60% for the central data. Even though our rates agree fairly well with the 
remaining data away from the p peak, there is still room for enhancement in 
the low mass region, 0.4 < M (GeV) < 0.6. 

We should finally mention what happens when the non-perturbative QGP 
rate is used instead of the perturbative result. Similar to the perturbative 
QGP results in Fig. 5.1.7 the non-perturbative plasma rate is about a factor 
or two larger in the low mass region. Even though this does help to explain 
some of the excess in the low mass region, especially in the more central data, 
the space time volume of the plasma phase is too small to have a large effect. 

5.1.8 Conclusions 

Using a parameterization of the results given by a hydrodynamic model of the 
collision region at the CERN SPS collider, the NA60 dimuon spectrum was re- 
produced using a pure thermal model assuming that there exists a sQGP phase 
above Tc with an interacting hadronic phase persisting until freeze-out. The 
dimuon spectrum from the sQGP phase originates primarily from qq annihi- 
lation with non-perturbative effects due to non-vanishing gluon condensates. 
After hadronization it is assumed that there remains a dilute hadronic gas 
in which the dimuon rates can be constrained entirely from chiral symmetry 
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Figure 5.5: NA60's excess dimuon data compared to our thermal yields which 
include contributions from either the perturbative or non-perturbative QGP 
phase, the hadronic phase and the DD contribution. Shown for all four cen- 
trality windows. 
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arguments and experimental data. The combination of these two rate equa- 
tions, after being folded over the space-time evolution of the collision region, 
are able to explain most of the excess dimuon data, especially in the more pe- 
ripheral collisions where our assumptions about diluteness hold. In the more 
central data, where the assumption of diluteness may breakdown, it is neces- 
sary to investigate how higher order terms in the virial expansion modify the 
spectrum. 

5.2 Di-electrons at RHIC 

In this section we now discuss the implications for the recent measurements 
by the PHEINIX collaboration at RHIC [118, 119]. Being a continuation of 
the work from the previous section (see also [13]) we do not go through all the 
details of the model for the hydrodynamic evolution. Table 5.2 summarizes 
the relevant parameters of the hydrodynamic simulation used in this section. 



Parameter 


Value 


EOS 


Lattice Motivated 


T c 


190 MeV 


T/.o. 


120 MeV 




0.2 fm/c 


impact parameter (b) 


f m 



Table 5.2: Parameters used in the hydrodynamic simulation of Au-Au colli- 
sions. 

There are some modifications that should be addressed. First, we now start 
the hydrodynamic evolution at tq = 0.2 fm/c in order to account for some of 
the pre-equilibrium dilepton production in the QGP stage of the evolution. 
This will increase the yields in the high mass region above the <p. It was 
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Figure 5.6: Summary of the equations of states used in this work. The lattice 
motivated EOS was taken from [120] and is used in the analysis of the RHIC 
data in this section. The BM EOS consists of a first order phase transition 
having a variable latent heat (shown for LH=1.2 GeV/fm 3 and 0.3 GeV/fm 3 .) 

already shown in the prior section (see also [10]) that non-equilibrium effects 
modify the dilepton spectrum but not the invariant mass spectrum. A 
second modification is the use of a more realistic equation of state motivated 
(taken from [120]) by recent lattice measurements. In figure 5.6 we show the 
squared speed of sound versus energy density for four different equation of 
states. The solid black curve is the lattice motivated EOS. The other three 
curves are bag model (BM) equation of states (i.e. first order phase transition) 
with variable latent heat. We show the BM EOS having a latent heat of 1.2 
GeV/fm 3 and 0.3 GeV/fm 3 . Also shown is the EOS having latent heat 1.2 
GeV/fm 3 with a fixed baryon number to entropy density (ub/s = 42). 

We now discuss the results of the space-time integration. In fig. 5. 7 we 
show the differential p± yields in four different mass windows. The thin lines 
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Figure 5.7: Differential spectra in different mass windows. The thin lines 
include the zeroth and first order, W£, contribution in pion density while the 
thick lines also included the second order, W£ n , contribution. 

show the zeroth and first order contribution, W n , and the solid lines include 
the additional two pion contribution, W n7T . We find a large enhancement from 
the two pion contribution at low p± which can be described as arising from 
Bremsstrahlung type processes. 

Figure 5.8 shows the invariant mass spectrum after going through the 
PHENIX acceptance. Shown separately are the contributions form the QGP 
and hadronic phases. Two scenarios are considered. The first is a hadronic 
phase with zero pion chemical potential. The second is the same hadronic evo- 
lution with a constant pion chemical potential of /i n = 50 MeV enhancing the 
rates at low mass throughout the entire evolution. The solid black line shows 
the sum of the contributions from the cocktail, QGP and hadronic phase at 

= 50 MeV. We find that after the acceptance cuts the one and two pion fi- 
nal state contributions are not able to explain the low mass enhancement seen 
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M (GeV) 

Figure 5.8: The data points show the measured di-electron spectrum from 
PHENIX. The dotted line shows the hadronic cocktail provided by PHENIX 
including the charm contribution. The solid curve labeled 'SUM' includes the 
cocktail, QGP and hadronic gas at p n = 50 MeV. 

in the data. In order to understand these cuts further we show in figure 5.9 the 
hadronic rates before and after the acceptance cuts. The yields in this figure 
are normalized at the p peak. It is clearly seen that the acceptance at PHENIX 
reduces the rates by about a factor of three or more for M < 0.5 GeV. To a 
good approximation, it turns out that the does not contribute at all after 
the acceptance cuts. Even though we find a large enhancement at low p± as 
seen in fig. 5.7, the sweet spot for the PHENIX detector is for dilepton pairs 
having m_|_ > 0.4 GeV while our 2ir enhancement feeds in below this kinematic 
region. In order to further the comparison with other works we also show 
the yields from 2ir annihilation proceeding through an intermediate p meson 
with a medium modified propagator. The calculation of the in-medium width 
is discussed in section 2.4.5. In figure 5.10 we show the hadronic rates from 
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Figure 5.9: Hadronic dilepton yields before and after the PHENIX acceptance. 
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Figure 5.10: PHENIX di-electron results compared to thermal 2tt annihilation 
with a collisionally broadened p. 

only 2tt annihilation with collisional broadening of the p propagator. 

Our results are consistent with the microscopic calculation of [121] when us- 
ing a collisionally modified rho showing that we at least have good constraints 
on the hydrodynamic evolution. Of course, we cannot explain any of the low 
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mass emission since there is no thermal emission below the 2ir threshold in 
this model. 

In summary, thermal hadronic emission helps explain the excess above the 
cocktail near the p mass. However, the low mass (M < 0.6 GeV) dilepton 
excess remains elusive. 

5.3 Role of Viscosity 

In this section we briefly discuss the effect that shear viscosity has on the dilep- 
ton spectrum. The motivation for looking at the viscous correction is twofold. 
First, there is the recently conjectured universal lower bound on the shear 
viscosity to entropy ratio, 77/s = 1/4-7T ~ 0.08 [48]. Therefore, one should see 
how viscosity modifies the current results which rely on kinetic equilibrium. 
Second, the empirical hadronic data on elliptic flow and pj_ spectra seem to 
support a small but non-zero shear viscosity. Since many of the simulations 
of the dilepton data at RHIC and CERN rely on kinetic equilibrium it will 
be interesting to see how off-equilibrium corrections modify the result. In- 
cluding viscosity will also help make contact with microscopic kinetic theory 
calculations. 

We will look at viscosity in the context of In-In collisions at NA60. Recently 
NA60 has measured the inverse slope (T e //) of the p± spectra as a function 
of invariant mass. Many of the theoretical calculations, including ours, were 
unable to account for the large rise in T e f / as a function of mass up to the about 
1 GeV. After 1 GeV there is a sudden drop in the effective temperature which 
is claimed by the NA60 collaboration to have origins from early emission before 
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the build up of collective flow thus resulting in its lower effective temperature. 

First let us look at our previous calculation (chapter 5.1) and see the effect 
of the changing the latent heat of the phase transition. Our previous calcula- 
tion used a latent heat of 1.2 GeV/fm 3 . Figure 5.3 compares the results for a 
latent heat of 1.2 GeV/fm 3 and a much smaller latent heat of 0.3 GeV/fm 3 . In 
terms of the mass spectra we find that the larger latent heat allows for more 
dilepton emission. This is simply because the space-time volume of the phase 
transition region is much larger. However, a large difference is seen in the 
T e ff results. We find that LH=0.3 leads to much higher temperature which is 
consistent with the data. There are two reasons for this large increase. First, 
a smaller latent heat is effectively a much harder EOS (see fig. 5.6). Second, 
differences in the space-time evolution lead to a larger fraction of the yield 
coming from the freeze-out contribution versus the hadronic contribution in 
the LH=0.3 EOS. 

Even though LH=0.3 does a decent job in explaining the T e ff data it will 
be hard to also explain hadronic observables in the same framework. It was 
shown in [35] that a latent heat of 0.8 GeV/fm 3 does the best job of fitting 
hadronic observables at the SPS. Second, this equation of state is far from 
consistent with the lattice motivated EOS of fig. 5.6. It turns out that the 
LH=0.8 and the lattice EOS lead to very similar effective temperatures in the 
dilepton spectra. 

With the above in mind we now examine how viscous corrections modify 
the dilepton spectra. We start by summarizing how shear viscosity modifies 
the particle / dilepton spectra by separating its role into two main effects. 
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Figure 5.11: The above figure shows the effect of changing the latent heat on 
the invariant mass spectra (top) and effective temperature (bottom) for central 
In-In collisions. Both hydrodynamic evolutions use T c = 165 MeV, To = 0.6 
fm/c and T/ . = 130 MeV. Also shown for comparison is the result using the 
Lattice EOS. 
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• Shear viscosity modifies the ideal equations of motion. For central colli- 
sion its main effect is to increase the transverse flow of the medium. This 
effect integrates over time and the largest change in the transverse flow 
profiles are seen at the latest times. The higher flow velocity will mani- 
fest itself by producing harder p±_ spectrum of the produced particles at 
freezeout. 

• Shear viscosity modifies the particle distribution function. The vis- 
cous correction to the distribution function schematically goes as 1/Tr. 
Therefore the viscous correction starts off large and becomes smaller as 
the system evolves in proper time. For realistic freezeout out surfaces 
it is found that the px spectrum hardens when including the viscous 
correction to the distribution function. 

One can now easily see how viscosity modifies the effective temperature. At 
early times the viscous correction to the distribution function causes an in- 
crease in T e ff while at later times it is instead the larger transverse flow which 
increases T e ff. This first effect will be most relevant during the QGP phase 
and the later effect most relevant during the hadronic phase. 

Let us now look at a quantitative example. We consider the hydrodynamic 
evolution of central In-In collisions starting at To = 1 fm/c. The QGP phase 
has a shear viscosity to entropy ratio of r]/s = 0.08 and we consider a hadronic 
phase both with (rj/s = 0.75) and without (rj/s = 10 -6 ) shear viscosity present. 
We also consider two scenarios for the dilepton production in the hadronic 
phase. In one case we consider the hadronic reaction nit — > p — > 7* — > e + e~ 
where the width of the intermediate p meson is collisionally broadened. In this 
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case one can include the viscous corrections to the rate in the same manner 
as was done for the qq case in chapter 4. We also do the calculation using the 
rates derived from the chiral reduction formalism in Chapter 2. However, as 
these rates relied on the use of kinetic equilibrium it is not trivial to include 
the viscous correction to the particle distribution function. Not including the 
viscous correction to the particle species is of course inconsistent, however, 
it is found that in the hadronic phase the largest viscous effects come from 
changes in the flow profile. Regardless, one can think of this result as showing 
the contribution from flow effects alone in the hadronic phase. 

In fig. 5.3 we show the effective temperature as a function of invariant mass. 
We still use a LH of 0.3 GeV/fm 3 but modify the freezeout temperature such 
that the maximum temperature agrees with the data at the rho mass. By com- 
paring with the ideal curves we can see how viscosity modifies the p± slopes in 
the two hadronic production scenarios. It is seen that even a modest viscosity 
in the QGP phase brings about a large increase in the temperature. Based 
on phenomenological cross sections [122] the shear viscosity in the hadronic 
phase is much larger. The increase in effective temperature is more modest 
since the gradients are smaller at later proper times. It is important to note 
that for M > 1.0 GeV the hadronic viscosity does not modify the tempera- 
ture since this region is dominated by qq annihilation in our model. At least 
in this model it appears that rj/s = 0.08 already saturates the experimental 
data. Any larger and the resultant temperature would be above the data for 
M > 1.0 GeV. 

Now we look at how the initial thermalization time effects the above result. 
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Figure 5.12: Effective temperature as a function of mass when the hadronic 
phase is constrained by the chiral reduction formula (top) or by thermal pro- 
duction through a medium modified rho (bottom). In each label the number 
in parenthesis is the viscosity in the hadronic phase. The upper curve in 
both plots is for r]/s qgp = 0.08 and rj/shad = 0.75. The middle curves are 
v/ s qgp = 0-08 an d v/ s had = 0. The lower lying curves are the results of the 
ideal simulation. All the hydrodynamic evolutions use T c = 160 MeV, tq = 1.0 
fm/c, T f . . = 135 MeV and a LH=0.3 GeV/fm 3 . 
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In order for the temperature to be consistent at the rho mass we have to in- 
crease the freezeout out temperature and LH of the model. Since changes in the 
LH and freezeout temperature don't effect the qgp temperature too strongly 
this wont effect the discussion to follow. The hydrodynamic evolution is now 
started at r = 0.5 fm/c. Focusing on M > 1 GeV the effective temperature 
of the ideal QGP phase is consistent with the data. However, when including 
even the minimal viscosity ofrj/s = 0.08 the temperature above M=l GeV is 
no longer consistent with the data. If we assume that the QGP phase must 
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Figure 5.13: Effective temperature as a function of mass when the hadronic 
phase is constrained by the chiral reduction formula. In each label the number 
in parenthesis is the viscosity in the hadronic phase. The upper curve is for 
v/ s qgp = 0-08 an d v/ s had = 0.75. The lower lying curve is the result of the 
ideal simulation. The parameters of the hydrodynamic evolution are T c = 160 
MeV, r = 0.5 fm/c, T f . . = 140 MeV and a LH=0.5 GeV/fm 3 . 

have a minimum viscosity of t]/s = 0.08 then we can conclude that within our 
model a thermalization time of Tq = 0.5 fm/c is inconsistent with the data 
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while a thermalization time closer to To = 1 fm/c is supported by the data. 

Finally, it should be discussed how the viscosity and thermalization time 
change the mass spectra. It was already shown in chapter 4 that the viscous 
correction to the distribution function leaves the mass spectra unchanged. 
However, changes in the flow profiles could possible change the spectra. Fig- 
ure 5.3 shows a summary of the mass spectra comparing a few different evolu- 
tion models used in the prior examples. To first order the mass spectra remain 
unchanged. 
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Figure 5.14: Summary plot of the invariant mass spectra compared to the 
NA60 data. The first four curves show the ideal and viscous {f]/s qgp = 
0.08,r]/ Shad = 0.75) using both the chiral reduction formula (xR-F) and col- 
lisional broadening (CB) with the parameter set described in fig. 5.3. The 
last two curves show the viscous and ideal curves using the parameter set of 
fig. 5.3. 
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Chapter 6 



Conclusions and Outlook 



Dilepton emission rates from a hadronic gas at finite temperature are con- 
strained by broken chiral symmetry in a density expansion. The rates are 
expressed in terms of measured e + e~ annihilation and r decay data. We focus 
on reactions having one pion and two pions in the final state. This brings 
about a mixing of the vector and axial spectral densities and results in a large 
enhancement of the rates below the 27r threshold. For the one pion state (corre- 
sponding to processes of the type / — * Fir+e + e~) this enhancement dominates 
at intermediate momentum (q > 0.5 GeV) and low mass (M < 2m n ). The 
two pion final state (/ — > Firn + e + e~) dominates at low momentum (q < 0.5 
GeV) and low mass. 

This result has direct consequences for dilepton phenomenology. Both 
NA60 and PHENIX see an increase in yield at low p±_ and low mass and 
we have shown that the 2ir final state contribution significantly enhances the 
yields in this kinematic region. Even though the enhancement can be thought 
of as arising due to Bremsstrahlung type processes it is a direct consequence 
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of the way chiral symmetry is broken. 

The second part of this work simulates a viscous relativistic hydrodynamic 
model of heavy ion collisions. We solve a second order theory that relaxes 
on short time scales to the Navier-Stokes stress tensor. Then the viscous 
correction to hadronic p± spectra and differential elliptic flow is calculated 
taking into account off-equilibrium corrections to the distribution function. 
In comparison to the ideal hydrodynamic results we find that shear viscosity 
decreases the differential elliptic flow and hardens the single particle spectra 
at larger transverse momentum. Both of these results are welcomed by the 
data. 

Next, the shear viscous correction to dilepton production from qq anni- 
hilation is computed. Starting from kinetic theory the quark's distribution 
function is corrected with its off-equilibrium counterpart. The resulting rates 
are integrated over the space-time evolution of the QGP phase using the above 
viscous hydrodynamic simulation. While the changes to the invariant mass 
spectra are small the viscous correction to the transverse momentum spectra 
becomes larger at higher invariant mass and/or p±. 

Given the success of ideal hydrodynamic models in their explanation of 
hadronic spectra a consistent description of electromagnetic observables should 
use the same space-time evolution. Therefore, in the last part of this thesis 
the dilepton rates are integrated over the space-time evolution of the hydro- 
dynamic evolution model. This is a large improvement over other models that 
employ parameterizations or phenomenologically motivated blast-wave type 
fits for the evolution. The results are then compared to the measured dilepton 
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spectra at NA60 and PHENIX. 

At NA60 the chiral reduction approach yields fairly good agreement with 
the invariant mass data. The agreement is better for more peripheral collisions 
with some discrepancy in the more central bins suggesting that the system 
might not be dilute enough for the density expansion. Part of the suppression 
at the p pole and the enhancement at lower mass which is seen in the data is 
a direct consequence of broken chiral symmetry. 

In regards to NA60's T e ff data, we have shown that the measured effective 
temperature can be obtained by either decreasing the latent heat, including 
viscosity in the QGP/Hadronic phases or some combination there of. If the 
mass region above 1 GeV is indeed dominated by QGP emission we have shown 
how the thermalization time and r]/s of the QGP can be extracted from the 
experimental p±_ data. 

Next we investigated the mass spectrum measured by PHENIX. Using 
the same comprehensive emission rates that worked on the percent level for 
NA60 the low mass enhancement seen at RHIC cannot be reconciled. This 
discrepancy can be traced back to the PHENIX acceptance which rejects all 
di-electron pairs having m± ^0.4 GeV. 

The future prospect for this work is to extract precision data, such as 
the EOS or transport coefficients, from the experimental heavy ion data. By 
fitting a combination of hadronic observables (p± spectra, elliptic flow) and 
electromagnetic observables (dilepton T e ff and mass spectra) within the same 
viscous hydrodynamic framework one could potentially extract quantities such 
as the thermalization time, latent heat, and shear viscosity to entropy ratio. 
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